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ABSTRACT 

A  threshold  detection  system  based  on  the  Neyman-Pearson 
criterion  is  derived  for  an  infrared  nutating  optical  syste 
Detection  is  optimum  for  small  signal-to-noise  ratios  and  a 
Gaussian  uncertainty  in  the  pointing  error  of  the  optical 
system.   The  signal  spectrum  and  background  power  spectral 
density  for  a  rectangular  space  filter  are  computed  numeri- 
cally and  used  to  specify  the  matched  filter  for  the  thresh 
old  detection  system  and  evaluate  its  performance. 
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I.   INTRODUCTION 

An  infrared  (IR)  detection  system  determines  the 
presence  or  absence  of  a  target  in  a  particular  area  of 
space  by  processing  the  received  IR  energy  from  that  area. 
The  target  must  be  found  amid  ever-present  background 
radiance  and  any  internal  noise  generated  by  the  detection 
system.   In  this  thesis,  one  type  of  system,  called  a 
nutating  optical  system,  is  considered. 

Infrared  detection  is  a  combination  of  spatial  frequency 
processing  and  temporal  processing.   Spatial  frequency  process 
ing  is  a  means  of  achieving  a  modulated  output  from  the  IR 
radiance  passing  into  the  optical  system  so  that  information 
may  be  easily  extracted  [Ref.  14].   To  extract  this  informa- 
tion, some  form  of  temporal  processing  is  used.   In  a 
nutating  optical  system,  the  optical  axis  rotates  about  an 
axis  perpendicular  to  the  image  plane  in  such  a  way  that  a 
point  image  traces  out  a  circle  in  the  image  plane.   The 
spatial  frequency  filter  is  a  piece  of  IR  sensitive  material 
in  the  image  plane  which  produces  a  voltage  output  that  is 
amplified  and  fed  into  a  temporal  filter  for  processing. 
A  nutating  optical  system  is  frequently  used  in  a  missile 
seeker  because  it  produces  a  signal  when  there  is  no 
tracking  error  and  the  rotating  optics  can  be  made  into  a 
gyro  that  provides  the  inertial  reference. 

This  work  is  a  synthesis  and  extension  of  two  recently 
published  articles  in  the  literature.   The  first,  by 


Harger  [Ref .  3] ,  formulates  the  detection  theory  of  a 
known  signal  in  background  and  white  noise.   Since  in  IR 
target  detection  the  signal  is  never  known  completely, 
Harger's  work  is  extended  to  the  case  where  the  target 
amplitude  and  position  are  unknown.   However,  a  known 
target  shape  is  assumed  as  well  as  a  priori  statistics  of 
amplitude  and  position.   Since  detection  is  most  difficult 
when  the  amplitude  is  small,  a  test  (called  the  threshold 
detector)  is  derived  that  is  optimum  for  weak  signals. 

To  specify  the  form  of  the  threshold  detector  and  to 
calculate  its  performance,  an  eigenvalue  integral  equation 
must  be  solved.   Ordinarily,  this  is  a  difficult  undertaking 
but  the  solution  becomes  trivial  if  the  covariance  function 
of  the  background  process  is  periodic.   This  is  exactly 
the  condition  realized  for  the  background  process  out  of 
a  nutating  detector  provided  that  the  phase,  which  is 
physically  unimportant,  is  averaged  out. 

Samuelsson's  results  [Ref.  11]  specify  the  form  of  the 
signal  and  background  coefficients  in  terms  of  the  system 
parameters  and  signal  and  background  radiance.   These 
derivations  were  first  published  in  a  Swedish  internal 
report  [Ref.  10],  then  re-derived  in  English  by  Yoshitani 
[Ref.  6].   The  latter's  derivations  are  included  in  the 
appendices . 

Samuelsson's  equations  have  been  applied  to  a  rectangu- 
lar IR  detector  assuming  a  Wiener  spectrum  for  the  random 
background  and  the  coefficients  have  been  calculated. 


Using  these  coefficients,  the  optimum  filter  has  been 
determined  and  the  probability  of  detection  computed  as 
a  function  of  signal-to-noise  ratios  and  background-to- 
noise  ratios. 


II.  SPATIAL  PROCESSING 

The  spatial  frequency  filtering  portion  of  the  detection 
system  is  comprised  of  a  lens  system  and  an  IR  sensitive 
detector  located  in  the  image  plane,   the  focal  plane 
of  the  optics.  To  adequately  derive  information 

from  the  IR  radiation,  the  lens  system  and  detector  must 
be  of  a  nature  that  causes  the  voltage  output  to  be  modulated 
in  a  way  that  information  about  a  target  within  the  field 
of  view  can  be  processed.   In  a  practical  nutating  optical 
system,  the  lens  system  moves  and  the  detector  is  stationary. 
However,  for  conceptual  purposes,  this  is  equivalent  to  a 
stationary  lens  systen  and  nutating  detector. 

A.   TARGET  SPECTRUM 

The  object  to  be  detected,  the  target,  is  located 

in  the  object  plane  a  distance  R  from  the  lens  and  parallel 

to  the  image  plane  as  shown  in  Figure  1.1.   The  angle  from 

the  perpendicular  axis  connecting  the  planes  to  the  target 

coordinates  measured  in  the  x-direction  is  called  x  and  is 

measured  in  milliradians ;  the  angle  in  the  y-direction  is 

y.   The  target  is  considered  to  be  close  to  the  optical 

axis  and  at  a  sufficient  distance,  R,  such  that  small  angle 

approximations  may  be  made, 


x  =  tan  x 

(2.1) 

y  =  tan  y 
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Consider  a  target  with  a  radiance  distribution 
N(r)  (W/cm  -sr)  which  is  the  radiant  power  into  a  unit  solid 
angle  per  unit  projected  area  of  the  sources.   r  =r(x,y)  is 
a  vector  whose  origin  is  the  perpendicular  axis  connecting 
the  coordinates . 

Assuming  no  atmospheric  attenuation  (this  can  be 
included  in  N(r)  if  it  exists)  and  perfect  optics,  the  power 
distribution  in  the  image  plane  is 


Pi(r)  -  N(r)  •  AQ  (w/sr) 


(2.2) 


where  A~  is  the  effective  entrance  area  of  the  optics. 

Dividing  by  A~  ,  one  can  see  that  the  radiance 
distribution  in  the  image,  plane  is  numerically  equal  to  the 
radiance  distribution  in  the  object  plane. 

However,  because  the  optics  is  not  perfect,  a 
point  object  given  by 


Hp  6(r)  =  Hp  6(x)  5  (y)  (w)  , 


(2.  ) 


(where  Hp  (W/cm  )  is  the  irradiance  at  the  optics  received 
from  the  point  source)  is  imaged  as  a  power  distribution 


Hp  •  fo(r) 


-1 


(W/cni   Sr) 


(2.4) 


where  f.(r)  (sr   )  is  the  point  spread  function  of  the 
optics.   Irradiance  is  defined  as  the  radiant  flux  indicent 
on  a  surface  of  unit  area  [Ref.  5]. 

Therefore,  the  radiance  distribution,  N(r),  in  the 
object  plane  will  be  imaged  as 

N'(r)  =  fN(r-s)  fQ(s)  d2s  (W/cm2Sr)  (2.5) 


12 


where  d  s  =  dxdy  and  the  integration  is  over  the  entire 
image  plane. 

The  detector  is  a  space  filter  in  the  image  plane 
which  may  be  a  moving  retical  or  simply  an  aperture  across 
which  the  power  is  integrated. 

The  power  in  the  image  plane  incident  on  the  detec 
tor  is 


H  =  Jn1 (r)  T(r)  d2r      (W/cm2) 


(2.6) 


where  i(r)  =  x(x,y)  is  the  transmi ttance  of  the  detector. 

When  the  detector  moves  as  a  function  of  time,  the 
power  incident  on  the  detector  becomes  a  function  of  time, 

H(t)  =jN'(r)  T(r-p(t))  d2r      (W/cm2)         (2.7) 

where  p(t)  describes  the  movement  of  the  detector  in  the 
image  plane. 

A  nutating  detector  moves  circularly  in  the  image 
plane  but  its  orientation  remains  fixed,  i.e.,  the  coordinate 
system  (x',y*)  shown  in  Fig.  1.2  does  not  rotate.   The  nuta- 
tion radius  measured  with  respect  to  the  optical  axis  is 
given  by  p. 

Assuming  a  nutating  optics  and  no  relative  motion 
between  the  target  and  the  optics,  the  radiance  on  the 
detector  will  be  periodic  and  thus  can  be  expanded  in  a 
Fourier  series 


H(t)  =  Z  H 


jnoo  t 
J      o 


0  <  t  <  T 


(2.8) 


n 
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where 


1      fT        -jnw  t 
Hn  =  i  J   H(t)  e 
o 


dt 


Replacing  H(t)  by  equation  (2.7), 


-.  rl  -jnw  t 

Hn  =  y  J   dt  e         d  J  N'(r)   (r-p(t)) 


(2.9) 


The  motion  describing  function,  p(t),  for  circular 
nutation  is 


p(t)  =  p(coswot,  sin  coQt)  , 


(2.10) 


where    oo     =    2ir/T    (rad/sec)    is    the    radian    frequency    of  nutation 
and   T    is    the    period    for    one    nutation. 

Substituting    (2.10)    in    (2.9),    H      can  be    shown    (see 
Appendix   A)    to   be 


H      =    jn    fN'(k)    T*(k)    e"3'11^   J    (2ttP  \lc      +k    2)    d2k  (2.11) 

nJ~~  n  xy 
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Figure  1.2.   Circular  Nutation  in  the  Image  Plane 
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where  k  =  (k  ,  k  ) ,  two-dimensional  spacial  Fourier 
x   y 

transform  coordinates  [Ref .  8]  , 

N'(k)  =  transform  of  radiance  distribution  in  the 

image  plane, 
T*(k)  =  complex  conjugate  of  the  transform  of  the 

optical  transmittance  function  in  the  detector 
coordinate  system, 
<J)  =  tan"   ky/kx> 
p  =  nutation  radius, 

and  J  (z)  is  the  nth  order  Bessel  function  of  the  first 
n 

kind.   It  has  been  assumed  that  the  image  plane  is  very 
large  compared  to  the  radiance  distribution  from  a  target 
so  that  the  image  plane  may  be  considered  infinite  for 
mathematical  purposes. 

B.   CORRELATION  OF  BACKGROUND  NOISE 

The  background  power  incident  on  the  detector  is  from 
a  sample  background  scene  where  a  scene  is  a  two-dimensional 
random  process  characterized  by  the  radius  vector  r  =   (x,y) 
from  the  optical  axis. 

Let  N(r)  be  a  sample  scene  radiance  distribution  on  the 
object  plane  from  an  ensemble  of  scenes  which  has  been  given  a 
suitable  probability  structure.   The  background  power  incident 

on  the  detector  is 

B(t)  =  /  N'(r)  T(r-p(t))  d  r  .  (2.12) 

The  covariance   function  is  defined  by  taking  the  expected 
value 
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E{[B(t)}  -  E{B(t)}]  [B(u)  -E{B(u)}]  =  Rg(t,u)     (2.13) 

where  E  denotes  the  expectation  over  the  ensemble  of  scenes. 
Because  the  detector  is  nutating,  the  random  process,  B(t), 
will  not  be  stationary.   However,  if  the  detector  and 
background  do  not  move  with  respect  to  each  other,  then 
RR(t,,t2)  is  doubly  periodic: 

RB(t1  +  MT,  t2  +  NT)  =  RB(t1}t2)  (2.14) 

where  M  and  N  are  integers. 

It  can  be  shown  [Ref.  7]  that  the  form  of  the  covari- 
ance  of  a  doubly  period  process  is 

oo       jmu)  t, 
m,n 

Because  the  nutation  is  periodic,  over  the  ensemble  of 
scenes,  there  is  little  significance  in  the  starting  time  t. 
Therefore,  considering  t  as  a  random  variable  that  is 
uniform  over  the  nutation  period,  one  can  "average  it  out" 
of  the  covariance  function  as 

1   T 
Rb(t)  «  i  /   RB(t+T,t)  dt  .  (2.16) 

Substituting  (2.16)  in  (2.15), 

0°       jmu    i   T   j  (m-n)  w  t 
Rb(t)  -   Z   8mn  e    °   1  /   e  dt.      (2.17) 

m,n  0 


16 


Using  the  relation 


i   T  j  (u-v)oont 
6uv  =  =r   /  e       u   dt  ,  (2.18) 

1  0 


the  correlation  function  becomes 


RB(x)  =  E  3n  e    u   .  (2.19) 

n 


Equation  (2.19)  shows  that  B(t)  can  be  considered  a 
wide  sense  stationary  process  and  that  it  is  periodic  in 
the  mean  square  sense;  i.e., 

Rb(t)  =  RB(T+x)   for  any  x  (2.20) 

The  power  spectral  density  is  the  temporal  Fourier 
transform  of  Rr(t).   Taking  the  transform  one  has 

oo 

SR(w)  =  2tt  E  6   5(w-ncon).  (2.21) 

B         n   n       0 

The  coefficients,  8  ,  can  be  related  to  the  optical 

'      n'  r 

system  parameters   by 


B      =    f|T(k)|2    |Fn(k)|2   WR(k)    J    2    (27rp/k    2  +  k  2)    d2k    (2.22) 
n        J'      ~    '       '    o   ~    '         B   ~        n  x        y  ~ 

where  Fn(k)  is  the  transform  of  the  point  spread  function, 
and  Wg(k)  is  the  Wiener  spectrum,  the  transform  of  the  back- 
ground correlation  function.   This  derivation  was  first  made 
by  Samuelsson  [Ref.  10]  in  an  internal  Swedish  report  then 
again  by  Yoshitani  [Ref.  6].   Yoshanti's  derivation  has 
been  included  as  Appendix  B. 
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Equation  (2.19)  implies  that  the  background,  B(t),  can 
be  represented  by  a  Fourier  series  with  uncorrelated 
coefficients  [Ref.  7].   Thus, 


B(t)  =  Z  bn  e 
n 


jnojgt 


(2.23) 


where  the  convergence  is  in  the  mean  square  sense 

The  coefficients,  b  ,  satisfy 

'   n '       J 


E(b  )  = 

*•  n' 


E(B(t)) 
0 


n  =  0 
n  +    0 


(2.24) 


Moreover,  they  are  uncorrelated 


E(b  b  ) 
v  m  n 


n 


0 


m  =  n 


m  f   n 


(2.25) 


C.   WHITE  NOISE  COMPONENT 

In  an  actual  system  there  will  be  a  certain  amount  of 
internal  noise  generated  in  the  detector  and  preamplifier 
which  is  additive  to  any  background  noise  or  signal.   In 
the  frequencies  of  interest  the  primary  source  of  this 
noise  is  thermal  agitation.   This  type  of  noise,  called 
Johnson  noise,  is  assumed  to  be  white  and  Gaussian. 

Let  W(t)  be  Gaussian  white  noise  with  power  spectral 

N    2 
density  y  (V  /H  ).   The  power  spectrum  at  the  output  of  the 

preamp  with  bandwidth  B  is  shown  in  Figure  1.3. 
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Figure  1.3.   Power  Spectrum  of  White  Noise. 


The  RMS  voltage  measured  at  the  output  of  the  preamp  is 


V2rms  =  j    •  2B  =  NB  (V2)  . 


(2.26) 
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Ill .   TEMPORAL  PROCESSING 

The  output  voltage  from  the  detector  is  temporally  pro- 
cessed to  determine  the  presence  of  a  target.   The  type  of 
processor  is  determined  by  applying  statistical  detection 
theory  to  the  known  behavior  of  the  nutating  system  and  var- 
ious realistic  assumptions  about  the  system  and  the  target. 
An  optimum  processor  in  the  sense  of  small  signal  amplitudes, 
called  a  threshold  detector,  is  developed  in  this  section. 

Harger  [Ref.  3]  derived  the  detection  theory  for  the 
known  signal  case  in  the  presence  of  additive  background  and 
white  noise  components.   In  IR  detection,  the  target  ampli- 
tude and  position  are  usually  not  known  completely.   The 
threshold  detector  is  an  extension  of  Harger' s  work  which 
accounts  for  the  unknown  amplitude  and  position. 

The  form  of  the  optimal  detector  under  the  Neyman-Pearson 
criterion  will  be  derived  for  multiple  observations.   It  is 
assumed  that  the  image  plane  and  object  plane  remain  parallel  and 
stationary  with  respect  to  each  other  during  the  observations 

The  detection  of  the  target  under  the  Neyman-Pearson 
criterion  becomes  a  problem  of  hypothesis  testing.   The 
task  of  the  detector  is  one  of  choosing  between  two 
hypotheses,  Hn  that  only  noise  is  present  and  H,  that  in 
addition  to  the  noise  there  is  a  target  present.   The  design 
of  the  detector  is  one  that  permits  the  correct  choice  of 
hypothesis  H,  (a  detection)  with  maximum  probability  of 
detection,  Qd,  in  a  fixed  probability,  QpA,  of  choosing 
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H,  when  Hfi  is  true  (a  false  alarm) .   The  structure  of  the 
Neyman-Pearson  criterion  requires  forming  a  likelihood  ratio, 
A,n,  which  is  compared  to  a  known  threshold,  V  ,  specified 
for  a  given  false  alarm  probability.   If  the  likelihood 
ratio  exceeds  the  threshold,  hypothesis  H,  is  assumed  true 
and  if  not,  hypothesis  H~  is  assumed  true.   The  logic  process 
may  be  written 

«i 

A10  <   V  (3.1) 

H0 

where  the  threshold,  V  ,  is  derived  from  the  expression 

CO 

QFA  =  /   p(Z|H1)  dZ  (3.2) 

Vt 

and  Z  is  a  sufficient  statistic  of  the  received  data. 

A.   DETECTION  OF  A  KNOWN  SIGNAL 

The  data  for  processing  is  a  set  of  M  functions  of  time 
{Z  (t) ,  n=l,...,M}  where  each  function  is  the  output  of  the 
preamp  during  one  nutation  (0,T).   T  is  the  nutation  period. 
The  output  may  be  represented  as 

Zn(t)  =  Sn(t)  +  B(t)  +Nn(t);  0  <  t<  T,  n=l,...,M   (3.3) 

S  (t)  represents  the  output  due  to  a  target.   The  component, 
B(t),  represents  an  additive  background  and  is  a  sample 
function  of  a  random  field  (B) .   The  relative  size  of  the 
target  in  the  field  of  view  is  assumed  sufficiently  small 
that  the  background  is  additive  while  B(t)  is  assumed  to  be 
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the  sample  function  during  each  of  the  M  nutations  observed. 
N  (t)  represents  the  additive  noise  component  due  to  the 
internal  noise  of  the  detector-preamp  combination.   It  is  a 
sample  function  of  a  white  Gaussian  random  process  (N  )  of 
zero  mean  and  power  spectral  density  N/2(V  /H?). 

The  two  hypotheses  for  the  likelihood  ratio  may  be 
written  as 

Hx:  {Zn(t  )  =  Sn(t  )  +  B(t  )  +  Nn(t  )  ;  n=l,...,M} 

(3.4) 
H0:  {?n( t  )  =  ?(t  )  +  NR(t  )  ;  n=l,...,M 

Denote  the  likelihood  ratio  as  A,n({Z  })  where  the  like- 

1U    n 

lihood  ratio  is  defined  as 

P({Z„}|H,) 

The  likelihood  ratio  is  complicated  by  the  fact  Z   and 
Z   are  not  independent  because  of  a  common  noise  component, 
B(t).   To  simplify  the  derivation  Harger  introduced  an 
auxiliary  hypothesis 

H-:  {Z   =  N   ;  n=l, . . . ,M}  (3.6) 

2    ~n    n  '     '    '   •  K 

Then  the  likelihood  ratio,  A-in,  can  be  computed  using 
the  chain  rule  for  likelihood  ratios: 

p({Zn>|H1)    p({Zn}|H1)/p({Zn}|H2) 
A10  =  p((Zn}|H0)  =  p({Zn} |H0)/p({Zn> IH^) 

=  A12/Aor 
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If  the  sample  functions  for  the  background,  B,  were  non 
random,  Z.(t)  and  Z . (t) ,  i  i    j,  would  be  mutually  indepen- 
dent because  of  white  noise.   Accordingly,  if  B  were  fixed, 
one  can  write  the  conditional  likelihood  for  the  n " 
observation  A, 2  ({Z  }|B).   Moreover,  the  likelihood  ratio 
for  M  observations  would  be  the  product  of  the  likelihood 
ratios  for  each  single  observation, 


M 


,  n 


A12  »?n»lBJ  =  \    A12  C?„l!> 

n=l 


The  conditioning  is  then  removed  by  averaging  over  the 
ensemble  of  backgrounds 


Since  the  noise  is  white  Gaussian,  the  conditional 

likelihood  ratio  is  that  for  detecting  a  known  signal 

S  +B   in  white  Gaussian  noise  and  has  the  well  known  form 
n   n 


([Ref .  15] ,  page  253) 


n 


A12  ^{?n}l!)  =  eXp 


?   M   T 

K     \    '   ZnW^(t)  ♦  B(t))  dt 
n=l  0 


1   M   T  2 

i  I      f       (S  ft)  +  B(t))Z  dt 

N  n=l  0    n 


which  may  be  rewritten  as 


n 


A12((Zn}|B)  =  exp 


?   M    T 

|  l     I     zn(t)  sn(t)  dt 

n=l  0 


1   M    T   2 

i   I   /   S'(t)  dt 

N  n=l  0   n 


exp  a(b)     (3.14) 


23 


where 


o   M   T  m  T  ? 

a(b)  =  £  E   /   (Z  ft)  +  S  (t))  B(t)dt  -  £  /  BZ(t)dt.     (3.15) 
1 n=l  0    n       n  iN  0 

To  find  A10^zn^>  one  must  evaluate  E  [An   ({Z  }|B)] 
where 


?     M   T 
EB[An12({Zn}|B)]  =  exp  {^   E   /  Zn(t)  Sn(t)  dt 

n=l  (J 


(3.16) 


1   M    T    2 
-  £  E   /   Snz(t)dt}  .  ER  {exp  a(b)}. 

n=l  0 


Assuming  the  background  is  Gaussian  the  expectation  can  be 
evaluated  by  representing  the  background  in  a  Karhunen- 
Loeve  (K-L)  expansion 


B(t)  =  Z  bR  g>k(t)  ,  0  <  t  <  T  (  .17) 

k 


where  ^(t)  is  the  eigenfunction  corresponding  to  eigenvaJ  le 
\,     of  the  homogeneous  integral  equation 


AkMt}  =  f      RB(t'u)  Mu)du  '  °  -  t  -  T'       (3.18) 


Rr>(t,u)  is  the  covariance  function  of  background 

RB(t,u)  =  E{[B(t)  -  mB(t)][B(u)  -  mB(u)]}  (3.19) 


where  mR(t)  is  the  mean  of  the  background. 


24 


Then  Eg(expa(b))  can  be  calculated  in  a  straightforward 

manner.   After  the  lengthy  calculation  shown  in  Appendix  C 

In  ER  {exp  a(b)}  =  -  1  Z  In  (1  +  J  A,) 

b  i    k         NQ      k 

?   T    M  T    M 

+  g-  /  dt  Z  (Zn(t)  -  Sn(t))  /  du  Z   (Zm(u)-S  (u))q(t,u) 

0    n=l  0   m=l 


T  M 

f  f   dt  (mB(t)  -  J  Z   (Zn(t)-Sn(t))}  h(t)         (3.20) 
0  n=l 


where  q(t,u)  satisfies  the  integral  equation 

T 

2g  q(t,u)  +  /   RB(t,T)q(T,u)dT  =  Rg(t,u)    0  <  u,  t  <  T   (3.21) 

and  h(t)  satisfies 

T 
^  h(t)  +  /   RB(t,i)h(T)dx  =  mB(t)    0  _<_  t  <  T.  (3.22) 

q(t,u)  is  the  impulse  response  of  the  filter  that  gives 
the  minimum  mean  square  error  estimate  B  of  the  background 
when  the  mean  is  zero  [Ref.  15].   The  additive  white  noise 
is  reduced  in  amplitude  by  1/M.   Using  the  K-L  representa- 
tion, one  can  also  write 


q(t,u)  =  Z  ^L_L_  q    (t)  ^*(u)   o  <  u,  t  <  T      (3.23) 

k  -1  +  X 
2M    Ak 


Using  equations  (2.9),  (3.10),  (3.14)  and  (3.20)  and 
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observing    that 


A02«Zn})    =    A12(Un})|Sn   =    0  (3.24) 


it  can  be  shown  (see  Appendix  C)  that  the  likelihood  function 
for  hypotheses  H.  and  H„  is 

M   T 

A   ({Z  })  =  exp  {  Z   /   Z  (t)q  (t)dt  +  K» }  (3.25) 

lu    n  n=1  Q    n     n 

where 

qn(t)  =  -  [Sn(t)  -  i  /    Z   S£(u)q(t,u)du]  (3.25a) 

U    X/ "~  J» 

and 


=    1   M   T    ?  TT  M    M 

*       N   Z   /   S   (t)dt  +  *±   ff    Z    Z   S  (t)S  (u)q(t,u)dtdu 
n=l  0  00n=l  £=1 

?   M   T 
-  -   Z   /   S  (t)h(t)dt  l  >.25b) 

n=l  0 


The  logarithm  is  a  monotonic  function.   Thus  it  offer 
no  loss  of  information  but  greatly  simplifies  the  calcula- 
tions.  In  the  remaining  work 


A'l0  ({Zn})  =  ln  A10  ({Zn})  *  (3*26) 


Rewriting  equation  (5.25), 

A'10({Zn})  =  L  ({Zn})  +  K'  (3.27) 

where 
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M   T 

L<  Zn  3  =   E.  '  Zn(t)qn(t)dt  . 
n=l  0 


K'  does  not  depend  on  the  input  data  so  it  may  be  included 
in  the  threshold.   The  hypotheses  test  now  becomes 

Hi 

L(  Z   )  J   Threshold  .  (3.28) 

Ho 

It  is  also  possible,  as  Harger  has  shown,  to  derive 
tests  similar  to  (3.28)  without  making  the  Gaussian  back- 
ground assumption.   In  particular,  if  the  background  fluctu- 
ates many  times  over  (0,  T)  such  that  a(b)  of  (3.15)  is 
approximately  Gaussian,  then  the  test  is 

Hl 

M   T  x 

I      f   Z  (t)g  (t)dt  „    Threshold  (3.29) 

n-1  0  n    n      < 

H0 


where 


i  ?m   T  l   M 


Since  by  Mercer's  theorem 


RB(t,u)  =  Z    Ak$k(t)  ♦*(u),    0  <  t,  u  <  T  (3.30) 


on  comparing  q  ,  g   and  (3.21),  it  is  seen  that  the  above 
model  approaches  the  Gaussian  model  if 
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k 

The  above  condition,  which  implies  a  weak  background,  is 
sufficient  but  by  no  means  necessary.   The  two  models 
approximate  each  other  whenever  the  second  term  of  q   and 
g   are  nearly  identical.   The  point  is  that  while  the 
Gaussian  background  assumption  may  not  fit  the  reality, 
a  processor  based  on  the  assumption  is  probably  not  far 
from  being  the  optimum.   Also,  the  Gaussian  assumption  allows 
one  to  compute  the  processor  performance. 

B.   DETECTION  OF  SIGNAL  WITH  UNKNOWN  PARAMETERS 

The  likelihood  ratio  derived  in  section  A  is  applicable 
only  if  the  signal  is  known  completely.   Such  is  seldom  the 
case  in  IR  target  detection  and  in  general  the  signal  is  a 
function  of  some  unknown  parameters.   Typically,  the  ampli- 
tude, a, of  the  target  and  its  position  r   in  the  image  plan 
are  the  unknown  parameters,  and  the  form  of  the  signal  can 
be  rewritten  to  include  these  parameters: 

S   =  S  (t;a,  r  ) 
n    nv  '  '  oJ 

Furthermore,  by  assuming  that  only  point  targets  are  of 
interest  and  including  the  assumption  of  fixed  image  and 
object  planes  over  the  number  of  nutations  of  interest, 
the  shape  of  the  waveform  will  be  the  same  during  each 
nutation.   However,  the  amplitude  is  allowed  to  vary  between 
successive  nutations. 
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The  signal  waveform  may  now  be  written  as 

Sn(t*a'V  =  an  f(t^o)  ^*31^ 


where  a   is  the  signal  amDlitude  and  f(t,r  )  is  the  signal 
n  ox  qj  o 

shape  for  a  target  at  r  . 

Rewriting  hypothesis  H,  under  these  assumptions 
gives 

H  '  :  {Z  (t)  =  a  £(t,r  )  +  B(t)  +  N  (t) ;  n=l,...,M 
1     n  ^  ^     n^'o^     y    J  ny    J  '     '    ' 

0  <  t  <  T}  (3.32) 


The  likelihood  ratio  A,fi  becomes  the  conditional 
likelihood  ratio  A,„({Z  }|a,r  )  obtained  by  replacing 
S  (t)  by  a  f(t,r  )  everywhere  in  (3.25). 

Assuming  a  priori  statistical  knowledge  of  a  and  r  , 
and  their  independence,  the  likelihood  ratio  may  be 
averaged  over  the  respective  density  functions  to  produce 
an  unconditional  likelihood  ratio 


A10({Zn})  =  JKo({V^Vp(a)p(VdadV     (3'33) 


Instead  of  averaging  as  above,  one  could  estimate  a  and  r 
by  the  maximum  likelihood  procedure  and  substitute  the 
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estimated  values  a  and  r„  in  the  conditional  likelihood 
ratio  and  try  to  maximize  the  generalized  likelihood  ratio 

A,~({Z  }  a.r,).   However,  this  procedure  is  not  followed 
10   n  ~ '  0  r 

here  because  the  maximum  likelihood  estimates  are  difficult 

to  obtain. 

Since  the  main  interest  is  in  detection  0f  the  weak 

target,  the  small  signal  case  is  of  primary  interest. 

Expanding  A,~({Z  )|a,r  )  in  powers  of  signal  amplitudes 

M  a 

A10({VLa>V    =    1    +    \    an  llT  A10^{Zn}l!'r0^a=0 

n  L  n  ~        (3.34) 

1  M   M        32 

+  i   £    Z   a  a,  a   -    A,n({Z  >la»rn)l   n 

2  ni1   n  k  3a  3a,   10v   n  '~   (r'a=0 
n=l  k=l        n   k 


Integrating  (3.34)  as  in  (3.33)  term  by  term  and 
retaining  only  the  first  non-zero  term  that  depends  on  the 
data  provides  a  statistic  which  is  optimum  for  weak  signals 
and  is  called  a  threshold  detector  [Ref .  4]  . 

The  evaluation  is  made  easier  by  the  fact 


3^  A10(tZn>l?>Vlar0  ^  ln   A10»Zn}l?>Vla  =  0         (3-35) 


which  results  from 


A({Zn)|0,r0)  =  1 
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The  first  non-zero  term  becomes 


■+■ 


sir  ln  W{zn}l2>VI  _n-  nfn  zn(t)f(t,?0)dt 

K  3.~  U       U 


£  /  Z  (t)dt  /  q(t,u)f(u,r  )du 
in  0   n       0  u 


/  f(t,r  )h(t)dt 
0       u 


(3.36) 


Performing  the  integration  over  r„, 

?      M      T  T 

Alf)((Z  })  »  1  +  jj  Z   a  •/  Z  (t)  [f(t)  -  /  q(t,u)  f(u)du] 
iU    n         1N  n  =  l   n  0   n  0 


/  f(t)  h(t)  dt 
0 


(3.37) 


where  the  bar  denotes  the  averaged  quantities.   The  first 
and  last  terms  are  independent  of  the  data  and  may  be 
included  in  the  threshold  value.   Moreover,  if  the  ampli- 
tudes are  identically  distributed,  as  is  generally  the  case, 
the  threshold  detector  performs  the  test 

(3.38) 


W  -   Z   /  Z  ft)  q(t)  dt  I      W 
11=1  0   n  „    Z 

H0 


where  W   is  the  threshold  and 


q(t)  -  |  [f(t)  -  /  q(t,u)  f(u)  du]  . 
1N         0 


(3.39) 
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Thus  the  threshold  detector  can  be  mechanized  either  as  a 
correlator  or  matched  filter  followed  by  an  accumulator. 

The  filter  function  q(t)  may  be  represented  in  a  series 
expansion  in  terms  of  the  eigenfunctions  of  RR(t,u)  and  the 
coefficients  of  the  signal  expansion.   Expanding  f(t)  in 
terms  of  the  eigenfunctions : 


f(t)  =  E  fk4>k(t)  •  (3.40) 


Substituting  (3.40)  in  (3.39)  and  using  (3.23) 


q(t)  =  £   Z  -*-  *  (t)  .  (3.41) 

k  2M  +  Xk 


C.   EVALUATION  OF  PERFORMANCE 

The  Gaussian  assumption  was  made  in  deriving  the 
threshold  statistic;  therefore,  the  threshold  statistic 
itself  is  Gaussian.   This  means  that  the  probability  of 
detection,  Qd ,  and  probability  of  false  alarm,  QF . ,  are 
uniquely  determined  by  second-order  statistics,  e.g., 
the  means  and  variances  of  W  under  hypotheses  Hn  and  H..  . 

Recall  that  the  threshold  statistic  is 

M    T 
W  =   Z   /  Z  ft)  q(t)  dt.  (3.42) 

n=l  0   n 

The  means  and  variances  of  W  are: 
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M    T 
E(W|H0)  =   I   /   q(t)  E(B(t)  +  N  (t))-dt 
n=l  0 


M    T 
=  I      f    q(t)  m  (t.)  dt 
n=l  0       B 


M    T 
E(W|H1)  =  I      f   q(t)  E[a   f(t,rQ)  +  B(t)  +  \(t)]  dt 
n=l  0 


M    T 
=  E(W|Hn)  +   E   /  f(t)  q(t)  dt 


0 


n=l  0 


VAR  (W|H0)  =  VAR  (W|H1)  =  VAR 

T  T      T 

=  ^T  /   q2(t)  dt  +  M2  /  q(t)  /  q(u)  RR(t,u)  dudt 
^  0  0      0       D 

(3.43) 

Letting  p  (»1H.)  be  the  Gaussian  density  function  of 
W  under  H.  which  is  specified  by  the  equations  above,  then 


!fa  =  '    Pw(xlHo'dx 


(3.44) 


and 


/  p  (xlH  )dx 


(3.45) 


where  W   is  the  threshold  voltage. 

Equations  (3.43)  and  (3.44)  can  be  solved  to  yield 
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and 


QFA  -  *c« 


(3.46) 


Qd  -  *c(x-d) 


(3.47) 


where 


x  = 


Wt  -  E[W|HQ] 

/VAR 


d  is  the  equivalent  signal-to-noise  ratio 


d  = 


EOVlH^ 
/VAR 


M      T. 

E   a   /  £(t)  q(t)  dt 
n=l   n  0 


J' 


(3.48) 


—  /   q2(t)  dt  +  M2  /   q(t)  /   q(u)  RR(t,u)  dudt 
0  0        0 


and  4)_(g)  is  the  complementary  error  function  of  the  kind 


+dCg) 


1    -   -v2/2 

—  f      e  dv 


Tn      g 


(3.49) 


By  setting  a  desired  probability  of  false  alarm,  x 
may  be  determined  from  equation  (3.46).   The  probability 
of  detection  is  then  determined  by  solving  for  d  and  using 
equation  (3.47). 

By  substituting  (3.41)  and  (5.42),  and  using  (3.23), 
the  numerator  of  d  is 
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M      T 

Z   a   /   f(t)q(t)dt 
n=l   n  0 


.  -  ,7 

M  f  * 

v    JL  v    '   K1 


"  2M  r  n 

n-l    k  2M  +  Ak 


(3.50a) 


while  the  denominator  of  d  (see  Appendix  D)  is 


/VAR  = 


_N_ 
2M 


2  1 


Z 
k 


2M    Ak 


Therefore  the  effective  signal-to-noise  ratio  is 


(3.50b) 


d  = 


M 


M 

I 

n=l 


n 


2M 
N 


1  + 


2M 


X. 


N   kJ 


(3.51) 


Knowledge  of  d  along  with  (3.46)  and  (3.47)  completely 
specifies  the  performance  of  the  threshold  detector. 

The  equivalent  signal-to-noise  ratio  for  a  target  whose 
amplitude  and  position  are  known  is  easily  shown  to  be 


df  = 


1   M 

—  T 

Mn=l 


II 


2M 

N 


ih 


k 


N   kJ 


(3.52) 


By  letting  the  background  be  zero  it  will  now  be  shown 
that  d'  reduces  to  the  equivalent  signal-to-noise  ratio  for 
a  known  signal  in  white  Gaussian  noise  which  is  well  known 
[Ref.  4] . 

The  energy  dissipated  during  one  observation  in  a  resis 

tor  of  1  ohm  if  the  signal  a  f(t,r  )  is  the  voltage  across 

&     n    '  o  ° 

that  resistor  is 
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T   °° 


OK  £ 


a2    Z    If,  (r    )  I  2 
n    ,     '    k v    cr  ' 


(3.53) 


The   voltage    resulting    from  M   observations    then    is 


—        1      M        — 

/E      =   n      E      ^E 

M  n-1        n 


1      *f 
th     £      a 

M  n=l     n 


J  MV 


(3.54) 


Sub 


stituting    (3.54)    in    (3.52)    withxk   =    0, 


d1    =    /E 


2M 
N 


2E      /M 


N 


which  shows  the  familiar  result  that  detectabil ity  depends 
only  upon  the  signal  energy  and  the  signal-to-noise  ratio 
increases  as  the  square  root  of  the  number  of  observations 
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IV.   THRESHOLD  DETECTOR  FOR  NUTATING  SYSTEM 

In  an  infrared  detection  system  the  output  voltage  to 

be  processed  is  obtained  from  a  detector-preamplifier  com- 

2 
bination.   The  irradiance  on  the  detector  in  W/cm   is  con- 
verted to  volts  with  a  linear  scale  factor.   The  voltage 
output  can  be  written  as 

v(t,rQ)  =  KH(t,rQ)  (4.1) 

-> 
where  H(t,r~)  is  the  irradiance  on  the  detector  from  a  target 

at  position  rn ,  and  K  is  a  scale  factor  in  (V/W/cm  ).   Also 

the  covariance  function  may  be  written  as 

RvB(t)  =  E[KB(t+x)  KB(t)]  =  K2Rb(t),  (4.2) 

where  B(t)  is  the  irradiance  on  the  detector  due  to  a  partic- 
ular background  scene.   The  output  voltage  due  to  this  see;  j 
is 

vB  =  KB(t)  (4.  ) 

The  problem  now  is  to  relate  the  detector-preamplifier 
output  to  the  equations  derived  for  the  threshold  detector. 

It  was  shown  in  Part  III  that  an  optimum  processor  for 
weak  signals  was  a  threshold  detector  that  performed  a 
certain  linear  operation  on  the  data.   The  linear  filter  was 
specified  in  terms  of  the  eigenvalues  and  eigenfunct ions  of 
the  integral  equation  (3.18),  where  the  kernal  RR(t,u)  is 
the  covariance  function  for  the  background. 
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In  Part  II,  the  covariance  function  was  noted  to  be 
periodic.   Rewriting  (4.2)  in  terms  of  the  output  voltage 

?  °°     jkco  (t-u) 

RvR(t,u)  =  K   I  Bk  e    °  (4.4) 

k 

Substituting  (4.4)  into  integral  equation  (3.18),  the  solu- 
tion for  the  eigenvalues  and  eigenfunctions  is  trivial  with 

jkco  t 
♦k(t)  - 


/T 


and 


AR  =  K2BkT  (4.5) 

where  T  is  the  nutation  period. 

The  signal  coefficients  were  also  noted  to  be  periodic 
and  could  be  expanded  in  a  Fourier  series 

00        j  kco  t 
v(t)  =  K  Z  Hk(r  )e    °  .  (4.6) 

k   K   ° 

Recognizing  the  similarity  between  equations  (4.6)  and  (3.41), 
the  signal  coefficients  for  a  known  signal  which  is  expanded 
using  the  basis  set  above  become 

f,  (r  )  =  KH,(r  )  /T~  (4.7) 

For  the  threshold  detector,  the  signal  is  averaged  with 
respect  to  an  a  priori  density  function  for  r^.   Averaging 
can  be  done  with  the  coefficients  to  give 
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->■  ^  .-* 


f,  =  /  p(r  )f,  (r  )dr 
k     i  v  oJ    k v  oJ      o 


-»-  ..  «->■ 


=  K/T  /  F,  (r  )p(r  )dr 
k  ^  cr x   cr   o 


(4.8) 


=  K/T  F, 


Using  equation  (2.10),  the  coefficients  Fv  can  be  written 


as 


H, 


jn  (dr  p(r  )  /N'(k,r  )t*  GO  e"  ^  J   2-rrp  Jk2  +  k2  d2k   (4.9) 
J      J      o^ K   o     ^  '  oJ        ~        n     \  x  y   ~   ^ 


The  transform  of  the  radiance  distribution  from  a  target, 
N'(k,r  )  is  a  function  of  the  target  radiance,  R(k),  the 
position  of  the  target,  r  ,  and  the  point  spread  function  of 
the  optics,  F  (k) .   Rewriting  (4.9), 


H, 


=    jn     fd2kR(k)    t*(k)    J      2-rrp "Jk2+kM  e"J"n* 


/ 


dr      p  (r    ) e 
o    *  K    GJ 


j  2-rrk  ,r 


(4.10) 


If  the  density  function,  p(r  ),  is  assumed  Gaussian  with  a 
standard  deviation  p  ,  then  the  averaged  coefficients  become 


H 


2   2   2 

r  -27TPz(k-  +  kz) 

JR(k)e  J 


x*(k)e";)'n<t)J   2ttd  Jk2  +  k2 
K-J  n 


d2k.   (4.11) 


k   J      J"W  L  ^^    unp"M\  x  y 

The  Gaussian  assumption  on  r   is  reasonable  because  r 

r         o  o 

often  represents  a  pointing  error  from  some  designated  tar- 
get position.   Then  p   describes  the  accuracy  of  the  pointing 
mechanism.   Substituting  (4.8)  into  (3.42),  the  filter  func- 
tion for  the  threshold  detector  becomes 

/ 

'k 


q(t)  =  K  I    Hv  /T 

k 


i  kco  t 

J         0 


[1    +   Q   K^3kT 


(4.12) 


/T 
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Rewriting  the  equivalent  signal-to-noise  ratio  in  terms 

of  equations  (4.8)  and  (4.12)  one  has 

M 

.,  1      a   /T/M 
M    ,   n 
n=l 


d  = 


/N7T 


k    J 


1  + 


MTK26k 


"N/2" 


(4.13) 


It  will  be  shown  in  Part  V  that  3,  may  be  represented  as 


Bk  =  aEVQ] 


where  aR   is  the  amplitude  of  the  background  noise  correla- 

2 
tion  function,  9,      is  the  instantaneous  f ield-of -view  of  the 

detector  and  Q-,  is  an  integral  function  of  the  optical  and 

detector  parameters  and  correlation  length  of  the  background 

radiance . 


Define  the  following  terms: 


Signal-to-noise  ratio  (SNR)  = 


1      M      -      - 

M        ,      n 
n=l 

/N/2 


(4.14a) 


a^K/T 


Background-to-noise    ratio    (BNR)    = 


/N/2 


(4.  4b) 


Using  (4.14a)  and  (4.14b),  d  may  be  rewritten  as 

-it 
2       1 


d  =  SNR  /M 


I      H, 


1  +  MQ,PNR' 


(4.15) 


Noise  equivalent  irradiance  (NEI)  may  be  defined  at  the 
detector-preamplifier  output  to  be 
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K-NEI  = 


7? 


2B 


/NB 


w 


(4.16) 


where  B   is  the  detector-preamplifier  bandwidth, 
w 

The  SNR  and  BNR  may  now  be  written  in  terms  of  the  NET  as 


/M  SNR  = 


<a>/NTB 
NEI 


SNR'  /NTB" 


where  NTB  is  number- time -bandwidth  product  and 


/M  BNR  =  rrp^-  /NTB  =  BNR'  /NTB. 
NEI 


Expressing  d  in  terms  of  NEI  one  has 


d  =  SNR'  /NTB  £ 
k 


H 


k 


1  +  Qk  NTB  BNR' 


(4.17) 
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V.   SOLUTION  FOR  RECTANGULAR  DETECTOR 

The  equations  in  Part  II  for  the  detector-preamplifier 
output  voltage  and  the  correlation  function  of  the  background 
must  be  solved  to  obtain  the  filter  function  and  specify  the 
detection  system  performance. 

The  two  general  equations  obtained  in  Part  II  are 

Hn(rQ)  =  3nJR(k)e"J"'r°Fo(k)Tni)e-jn* 


J[2ttp  /k2+k2  |d2 

where    (r  )  is  the  n    Fourier  series  coefficient  for  the 
radiation  on- the  spatial  filter  and 

K   -   /l^)|2    lFo^|2   W    JnPPVkx+ky  Vk 

where  B   is  the  n    coefficient  of  the  background  correlation 
n  ° 

function  at  the  spatial  filter. 

A  rectangular  sTiaped  detector  was  chosen  for  the  calcu- 
lations because  it  exhibited  the  simplest  form  of  the 
equations.   For  circular  nutation  the  filter  is  represented 
graphically  in  Figure  5.1. 
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Figure  5.1.   Nutation  with  Rectangular  Detector. 


If  the  detector  function  is  given  by 


T(x',y') 


x  -w/2  <  x*  <  x  +w/2 

1    0      —  —   o 

y0-h/2  <  y*  <_   y0+h/2 

0   otherwise 


then  its  transform  is 


x(k  ,k  )  = 


-  i  2 tt  ( k  x  +k  v  ) 
J   ^  x  o  y  oJ 

Ti2k  k 
x  y 


sin(Trk  w)sin(irk  h) 
x        y 


so  that 


T(k) 


2  2 

sin  (-rrk  w)   sin  (irk  hj 

j\.  y 


7T  k 
X 


1,1 

TT  ky 


(5.1) 


The  point  spread  function  of  the  optics  is  assumed  to 
be  Gaussian: 
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I    1 2 

fo(r)    -*,      ^ 

2TTG 

Its    transform   is 

-2tt2q2    (k2+k2) 
F0QO    =   e  X      y        ,  (5.2) 

The  correlation  of  the  background  radiance  is  assumed  to 

be  of  the  form 

j.  f         \  2   -a  I  x  [  -  8  ly  I  ^^■7^ 

<l>B(x,y)  =  aR  e   '  '  e  p  ^   '  (5.3a) 

where  a    and  8    are  the  correlation  lengths  in  the  x-  and 
y-  directions,  respectively.   The  Wiener  spectrum  is  the 
Fourier  transform  of  <+>~: 

2      2a           26 
WR(k)  =  aR  -= =■  -* — *-  (5.3b) 

*  ~      C  a   +  (2-rrk  )Z  gZ  +  (2irk  )Z 

Substituting    (5.1),    (5.2)    and    (5.3b)    into    (2.10),    the 
radiance    from   a  point    target   becomes 

r  -i27r(k    r    +k    r    )       -  2tt2q2  (k2  +  k2) 
li    r~*   \         -n    /      J       *■   x   x     y   y  v   x     yJ 

H(r)=j/e  7/e  7 

n  *•    oJ         J     J 

i  2tt  (k   x   +k   y    )    sinuk  w   siniTk   h 
J       K   x   o      y7o^  x  y 

irk  irk 

x  y 


-jn   tan^Cky/k^     t    £_       r^^ 


J      2-^p    /k"+k"    dk   dk    .       (5.4) 
n\         v    x      y        x      y 


Letting 


U         =      /2  TT  G  k 
X  X 


Uy     =      /2TTOky 
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equation  (5.4)  becomes 


^      .n  r   - (u  +u  ) 

TT  J 


-j —  |(r  -x  )u  + (r  -y  )u  1 

o  L  x   o;  x  ^  y  ;o;  yj 


-jn  tan   (u  /u  ) 


y/uxJT   /2p   /  2A  2 

;  3      — -   /u  +u    du  du     (5.5) 

n   a  \l    x   y  I   x   y 

If  r   is  random,  the  radiance  from  a  target  may  be 
averaged  with  respect  to  the  density  function  p(r  )  as  shown 
in  Part  IV.   In  particular,  if  r   is  Gaussian  distributed, 
with  standard  deviation  p   the  averaged  coefficients  become 


ft.-JnJ- 


-2TT2Y2(k2  +  k2)      j27i(k   x   +k  y   ) 
1     v    x      yJ      J       ^   x   o      y7o^ 


sin-nk   w  sinuk  h        -in    tan      (k  r/k    ) 

x  y      e  J                  y     xJ 

Trk  uk 

x  y 


J      I  2irp    /k2+k2    |dk   dk  (5.( 

n  J    x     y  x     y 


where 


2  2  2 

Y      «   P0    +    o    . 


Letting 


U  '     =      /2~TTYk 

x  '    x 


Uy      =      /2TTYky 


equation    (5.6)    becomes 
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H 


n 


.j]    r  -(u2-u2)      j  —  (u'x   +u'y    )       -jn    tan-1  (ii'Vi/') 
2_    J  x      n        Jy^xo      y7  o  J  x      y 

2    /e  e  e 

IT       J 


J 


n         y 


/2P    rTpi 

;u    +  u        du 'du  ' 


x      y    I      x      y 


(5.7) 


And   finally   substituting    equations    (5.1),    (5.2)    and 
(5.3b)    into    equation    (2.22),    the    correlation    coefficients 
become 


6n    "      e 


2    2      2      2              2  2 

-4tt    a    (k    +k    )    sin    (ttV.   w)    sin    (irk   h) 
v    x      y7    __2L_  y 


(Trk    ) 
v      xy 


(irk    ) 

^     yJ 


2a 
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B   a2    +    (2irkv)2    B2    +    (2Trk) 
a  y 


Letting 


-  J2  [2ttP    /k2  +  k2  jdk   dk 
n  1         \J    x      y  I      x      y 


(5.8) 


u"   =    2-rrpk, 


Uy     =      2^P\ 


equation  (5.8)  becomes 


2n2     '    2  '2    "2    /Sin 

oDfi    a3o       r  -u    -u      / 


TT 


U"   \ 
X  \ 


sm 


u' 


U" 

/la      x 


/2"< 


v 


/2a     y 


2     '2 
(aa)    +u      (8a  Y'  +  u 
v      J         x    K  y 


1 j2   IP 

7    "2      n      a 


ti?    •  »?    \ 

u   +u       idu"du" 
x      y  x     y 


(5.9) 


where 


n2  2,  2 

fi      =  w  h    . 


(5.9) 
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The  integrals  in  (5.5),  (5.7)  and  (5.9)  were  difficult 
to  evaluate.   Although  the  components  of  the  integrand  are 
relatively  simple  functions,  the  argument  of  the  Besse] 
function  prohibited  separating  the  integral  into  two  one- 
dimensional  integrals.   No  closed  form  solutions  were  found 
so  the  integrals  were  evaluated  numerically.   The  numerical 
method  used  is  based  on  work  by  Pierce  [9]  who  applied 
Gaussian  quadrature  formulas  to  two  dimensional  integration 
by  integrating  over  a  planar  annulus  in  the  (x,y)  plane. 

Gaussian  quadrature  formulas  are  a  means  of  evaluating 
the  integral  by  summing  weighted  values  of  the  integrand  at 
specific  points.   For  the  one  dimensional  case, 

b  b 

/   g(x)dx  =  /   w(x)f(x)dx  (5.10) 

a  a 

where  g(x)  is  the  integrand  to  be  integrated  and  w(x)  is  a 
weight  function  for  which  the  specific  integration  was 
derived.  ,  For  well  behaved  functions  the  integral  may  be 
evaluated  as 

b   •         m 
/   g(x)dx  =   E   A.f(x.)  +  error  (5.11) 

a  i=l  1        1 

where  the  x.'s  are  the  m  zeros  of  the  m   polynomial 

m 

p  (x)  =  n  (x-x.) 

i  =  l 

of  the  set  of  polynomials  mutually  orthogonal  on  the  inter- 
val [a,b]  with  respect  to  the  weight  function  w ( x) .   The 
weights  are  given  by 
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Ai  =  /   w(x)  Li(x)dx  (5.13) 

a 

where  L^(x)  =  Pm  (x) / [  (x-x,)P '  (x . ) ]  is  the  Lagrange  inter- 
polation coefficient. 

For  a  weight  function  w(x)  and  interval  [a,b]  the  poly- 
nomial P  (x)  and  its  zeros  x.  and  weight  factors  A.  need  be 
m v  J  1        to  1 

computed  only  once.   For  a  number  of  weight  functions  and 
intervals  the  set  of  orthogonal  polynomials  is  known.   Stroud 
and  Secrest  [13]  give  x-'s  and  A.'s  for  a  variety  of  weight 
functions  and  internals.   The  degree  of  the  formula  deter- 
mines the  number  of  x. 's  and  A- 's.   A  2M-1  degree  formula 

11  b 

will  have  M  points  and  is  exact  for  polynomial  integrands  of 
degree  2M-1  or  less. 

Pierce  applied  Gaussian  quadrature  integration  to  two 
dimensions  where  the  solution  is  of  the  form 

//g(x,y)dxdy  -  E  E  D. ,g (x.. ,y.. ) 

i  j    J     J         J 

The  integration  is  over  an  annulus  in  the  x,y  plane  wi  h 
inner  radius  R  and  outer  radius  1.   Rewriting  the  integral 
in  polar  coordinates 

1  2tt 
/  /    rG(r,0)dedr  *  Z  I    D..G(r-,6.)  (5.14) 

R  o  i  j    J    J 

Pierce  showed  the  summation  could  be  rewritten  as 


4(m+l)  m+1 
E      E   A-B.G(r.  cos9.,r-  sine.)  . 


(5.15) 
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where  A^  and  R.  are  the  weights  for  the  radial  and  angular 
directions  respectively  and  4M+3  is  the  degree  of  accuracy. 
For  an  arbitrary  annulus  the  formula  may  be  obtained  by 
rewriting  equation  (5.14) 

2tt   r2 
I  =  /    /    f(r,0)  rdrde  (5.16) 

0   rx 

and  letting 

r 

P  -  ~ 

r2 

?   2tt   1 
I  =  rj;  /    /      pf(r2p,6)  dpde.  (5.17) 

o   r,  / 

l/r2 

The  approximate  solution  can  be  derived  as 


I  = 


4 (m+1)  m+1       /  r^ —2 n  \ 


where 


A.  =  27T/4(m+l) 

and 

2   2 

B.  =  w.  (ri;-r,  )/2,   w.  =  weight 

M+1  is  the  order  of  the  orthogonal  polynomial,  in  this  case 
the  Legendre  polynomial  on  the  interval  (0,1)  anda.  's  are 
the  zeros  of  this  polynomial.   This  type  of  formula  is  known 
as  a  spherical  product  Gauss  formula  [12]. 

The  annuii  were  picked  at  radii  corresponding  to  the 
zeros  of  J  (z)  because  J  (z)  is  the  fastest  damping  component 
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of  the  integrands.   One  additional  annulus  inside  the  first 

2 


zero  was  found  necessary  to  account  for  the  T   term  when 


a+x 
a  is  small  and  x  (and  y)  approach  zero.   The  distances 

between  annuli  outside  the  first  five  to  eight  zeros  were 

found  not  to  be  critical  and  thus  spaced  arbitrarily. 

The  Gaussian  Quadrature  formula  used  was  a  Gauss -Legendre 

24-point  formula  as  listed  in  [13]  thus  giving  a  degree  of 

accuracy  of  99.   The  integration  is  limited  to  only  the  first 

quadrant  as  mentioned  below.   Thus  for  a  24-point  formula, 

2 
(M+l)   or  576  points  per  annulus  were  used  to  evaluate  the 

integral . 

The  expression  for  3   is  easily  seen  to  be  an  even  func- 
tion both  radially  and  about  both  planer  axes,  therefore, 
only  integration  over  the  first  quadrant  was  necessary. 

With  a  little  more  difficulty,  F  (r  )    and  H   can  be  shown 

1  '   n  *•  o^       n 

(see  Appendix  E)  to  exhibit  the  same  property  provided  one 
expression  is  used  for  even  n  and  another  for  odd  n.   By 
calculating  only  over  the  first  quadrant  computer  time  was 
considerably  shortened. 

In  order  to  solve  for  coefficients  of  arbitrary  order, 
the  Bessel  functions  of  that  order  must  first  be  obtained. 
The  available  library  subroutines  were  found  inadequate  for 
this  use  for  two  reasons.   First,  these  routines  have  an 
upper  limit  of  100  on  the  order  of  the  Bessel  function  that 
can  be  computed  and  second,  the  routine  had  to  be  called  for 
each  individual  order. 
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The  recursive  relation 


J   (y) 

n  w  J 


2n 

y 


Vl^  -  Jn+2^ 


(5.19) 


provides  a  basis  for  generating  an  array  of  Bessel  functions 

of  various  orders.   However  for  accuracy  J  , , (v)  and  J   - (y) 

J      n  +  1  '  J  n+2  w  J 

must  be  known. 

The  most  accurate  method  to  numerically  calculate  Bessel 
functions  of  various  orders  and  arguments  was  found  to  be  a 
uniform  asymptotic  expansion  involving  Airy  functions  [1]: 

(A) 


t  r       i      4A 


*-       2  Z7, 

T71    J"„  ~~2l 


ii 


k=0   n' 


A»(n2/3A)   co   K(A) 

+  _± y     K 

n      k=0   n 


(5.20) 


where  A.  and  A!  are  two  types  of  Airy  functions  and 


2  A3/2    „ 

3  A     =  £n 


I  I  / 21/ 

1  +  Vl-z^  \/z 


JT7 


z  <  1 


aQ(A)  -  1 


b  (A)  =  -5/48  A2  +  — 
0  /A 


24(l-z2)3/2 


8(l-z2)'2] 


<  1 


av(A)  <<  a  fA) 


bR(A)  <<  bQ(A) 


k  >  1 


(5.21) 


and 


The  Airy  functions  are  calculated  using  the  equations 


A^x)  =  \  x~h   e"5f(-6) 


A!(x)  =  -  \  xh   e"6g(-6) 
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where 


and 


6  =  |  x3/2  (5.22) 


2/3  . 
x  =  n    A 


The  functions  f(-6)  and  g(-o)  are  tabulated  and  linear 
interpolation  is  used  to  determine  their  values.   A  is  guar- 
anteed to  always  be  positive  if  n  is  selected  to  be  larger 
than  the  argument. 

By  picking  an  order  much  higher  than  the  argument,  v, 
and  J   ,  ,  using  the  asymptotic  expansion,  the  recursive 
relationship  (5.19)  may  be  used  to  generate  an  array  of 
successive  orders  of  the  Bessel  function  down  to  and  includ- 
ing J  .   However,  error  builds  up  rapidly  using  this  method. 
A  second  relation 

1  =  JQ(y)  +  2J2(y)  +  2J4(y)  +  ...  (5.23) 

may  be  used  to  generate  a  normalizing  factor 

C  -  l/[j0(y)  +  2J2(y)  +  2J4(y)  +  .  .  .]  (5.  4) 

Multiplying  each  Bessel  function  value  generated  in  the 
recursive  relation  by  C  resulted  in  very  accurate  values 
being  obtained.   The  computer  program  used  to  solve  the 
equations  is  listed  as  Appendix  H. 

To  evaluate  the  numerical  results  a  check  was  used  which 
summed  the  coefficients  and  compared  this  summation  to  an 
analytic  expression  for  a  sum  that  could  easily  be  solved 
numerically . 
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For  a  particular  phase  angle  of  the  nutation  cycle,  co  t, 
the  radiance  function  for  a  target  is 


Hfni  t)  s  Z  H,  (r  )e 
K   o    J         ,   k  *•  o^ 


j  kco  t 
J      o 


(5.25) 


For  a  point  target  and  a  rectangular,  circularly  nutat- 
ing detector 

2  2  2 
-2tt  a  k   sinirwk    j  2irk  (pcosto  t  -  (r  -x  )) 

H(oj   t)    =   fe  x  r- — -  e  x  °  x     °     dkY 

v    o    J  Trk  x 

x 

2    2    2 
-2tt    a    k      sinnhk        j27rk    (psinto    t-    (r    -y    )) 

-fe  y   -_ X  e    y      °     y  °   dk 

uk  y 


(5.26) 


which  has  the  solution  (see  Appendix  F) 


H(o)  t)  =  t 

v  o  J         4 


E,  (r  ,x  .w)  -  F,  (r  ,x  -w) 


1 v  x'  o 


1  v  X*  o' 


-  E2(ry,y0,-h) 


EJr  ,y  h) 

2  v  y ''  o  ' 


where 


E1(r,z,a)  =  Erf 


r r- z -pcosw  t 

o_     a 


/2a 


2/2~a_ 


E2(r,z,a)  =  Erf 


r  r-z  -psmti)  t 

o      a 


/la       2/2a_ 


and  Erf  is  the  error  function 


7   x   _  2 
Erf(x)  =  —  /   e     du  . 
/?  0 


(5.27) 


By  evaluating  the  analytic  expression  at  the  phase  in- 
stant (w  t)  that  the  detector  crosses  the  target  and  comparing 
o 


55 


it  against  equation  (5.27)  for  an  appreciably  high  n,  a 
reasonable  check  was  made  on  the  accuracy  of  the  individual 
signal  coefficients. 

For  the  coefficients  averaged  over  r  ,  the  check  was 

°  o 

made  at  w  t  =  0  where 
o 


H(o)  -  I  H, 


(5.28) 


and  the  analytic  expression  is 


-  E3CVy0,h)] 


where 


E3(r,z,a)  =  F2  (r,z,a) 


a)  t=0 
o 


(5.29) 


The  rms  background  voltage  was  used  to  check  the  coef- 
ficients in  the  correlation  function.  The  mean  square  vol 
age  is 


eg  =  E{B2(t)}  =  RB(0)  =  I    6k  . 


(5.30) 


Appendix  B  shows  that  the  correlation  function  may  be  written 


as 


2  i.i  t 


O      ,Z- 


Rb(t)  =  /  x(k)r  W'(k)  J   4^pk  sin  -±-  d"k     (16B) 

Substituting  (5.1),  (5.2),  and  (5.36)  in  (16B)  and  evaluat- 
ing at  x  =  0 
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eR   =   o-p.2a   f- 


sin    (irk   w) 


-    2    2,  2 
4  it   a    k 
x 


Okx)2      a2    +     (2rrkx) 


T  dk 
2         x 


sm    (uk   h)         e  y 

aB2e  f     ,  ,    J        72 T^TT 7 


(irk    ) 


bz  +  (2TTky)z     y 


(5.31) 


It    can   be    shown    (see   Appendix    F)    that    an    analytic    expres 
sion    is 


e2    =    4aBa2R(w;a)R(h;3) 


where 


R(p,a)    = 


a 


1-Erfc 


I 


2y 
a 


2  a 
/tTij 


1  -  e 


l2aJ 


2    2 
a   a 


2a' 


e    M   Erfc    aa    +    ■—- 
2a 


-ay    r    r 
+    e  Erfc 


aa    -    J—      -    2Erfc(aa) 
2a 


and   Erfc    is    the    complimentary   error    function 


(5.32) 


?   °°    2 
Erfc(x)  =  —  /   e"u  du  . 

/F  x 


(5.33) 


The  computer  program  used  to  solve  equations  (5.6), 
(5.7)  and  (5.9)  and  calculate  the  infinite  summations  is 
listed  as  Appendix  H. 
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VI .   NUMERICAL  RESULTS 

Numerical  results  constitute  solving  for  the  signal  and 
background  coefficients  and  using  these  coefficients  to 
specify  the  form  of  the  filter  to  be  used  in  the  threshold 
detector  as  well  as  to  calculate  the  probability  of  detection 
as  a  function  of  signal-to-noise  ratio  (SNR) ,  background-to- 
noise  ratio  (BNR)  and  number  of  nutations  (M)  used  to  make 
the  decision.   The  probability  of  detection  is  determined 
for  the  averaged  equivalent  signal-to-noise  ratio  and  for 
specific  point  targets. 

Many  different  values  for  system  parameters  and  target 

locations  were  used  to  compute  the  coefficients.   A  typical 

set  of  parameters  is  : 

nutation  radius  =  15 

width  of  detector  =  30 

height  of  detector  =  1 

blur  circle  standard  deviation  =  .5 

background  correlation  length  x-direction  =  20 

background  correlation  length  y-direction  =  20 

position  of  detector  in  detector  coordinate  system 

x  direction  =  15 

y  direction  =  0 
target  coordinates  (if  desired) 

x  direction  =  0 

y  direction  =  0 
pointing  error  standard  deviation  =  7.5 

A  detector  utilizing  the  parameters  above  would  trace  out 

the  area  in  the  image  plane  shown  in  Figure  6.1. 


All  units  are  in  milliradians 
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t  =  T/2 


wrnm^i 


Figure  6.1.   Scanned  Area  of  Nutating  Detector. 

This  type  scan  was  chosen  because  it  crossed  the  center 
of  the  coordinate  system  which  would  be  the  center  of  the 
target  blur  circle  for  zero  pointing  error.   Also,  in  the 
case  of  a  random  pointing  error  in  the  absence  of  bias  ze  3 
mean  is  the  center  of  the  coordinates.   This  is  important 
because  the  threshold  detector  design  is  based  on  a  Gaussian 
pointing  error  with  zero  mean. 

A  sample  calculation  of  121  coefficients  using  the  para- 
meters above  is  listed  in  Appendix  G.   The  equations  shown 
in  Part  V  were  used  to  check  the  accuracy  of  the   coefficients 
Table  6.1  lists  typical  values  for  the  summation  of  121  coef- 
ficients for  various  sets  of  parameters  and  the  calculated 
summation  values. 
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TABLE  6.1 


Calculated  Sum     Actual  Sum 


Signal  (point  target)       .6826894  .6826993 

.6826895  .6826994 

.5438002  .5438087 

Signal  (averaged)  .2990664D-6  .2991232D-6 

.2328138D-4  .2328142D-4 

.2651595D-1  .2651595P-1 

Background  .6013477  .6012230 

.2347420  .2347420 

.8335525  .8319764 


Besides  the  general  shapes  of  the  background  and  signal 
spectra,  variations  in  the  spectra  due  to  detector  size, 
background  correlation  lengths  and  coordinates  of  the  detec- 
tor are  of  interest.   Curves  in  the  following  figures  show 
the  effects  of  these  variations. 

Figures  6.2  through  6.5  show  the  envelope  of  the  one 
side  of  the  two-sided  background  power  spectral  density 

00 

SB(w)  =  Z    6k6(oj-ku)o) 
k 

where  uj   is  normalized  to  2tt.   The  curves  of  Figures  6.2  and 
o 

6.3  illustrate  the  effect  of  varying  the  background  correla- 
tion lengths.  In  Figure  6.2  the  correlation  lengths  are  the 
same  in  both  x-  and  y-directions .  Practical  measurements, 
however,  indicate  that  correlation  lengths  in  the  horizontal 
(x-)  direction  may  be  longer  than  in  the  vertical  (y-)  direc 
tion.   The  curves  in  Figure  6.3  are  for  this  case. 
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The  curves  of  Figure  6.4  and  6.5  relate  the  spectrum  to 
changes  in  size  of  the  detector.   In  Figure  6.4,  only  the 
detector  with  BNR  equal  to  9  and  number  of  nutations  equal 
to  100  are  shown  in  Figures  6.11  through  6.13.   Each  figure 
shows  the  envelope  of  one  side  of  a  double  sided  spectrum. 
Two  cases  are  considered.   First,  threshold  detection  systems 
designed  for  different  Gaussian  pointing  errors  but  with  the 
same  size  detector  are  shown  in  Figures  6.11  and  6.12  and 
second,  detection  systems  designed  for  the  same  printing 
error  but  utilizing  different  size  detectors  are  shown  in 
Figures  6.12  and  6.13. 

The  curves  shown  in  Figures  6.14  through  6.16  show  the 
matched  filter  to  a  particular  point  target  located  within 
the  scan.   Notice  each  of  the  filters  exhibits  a  band  pass 
characteristic.   This  is  because  of  the  additional  noise 
component  (background) .   The  oscillations  observed  in  the 
point  target  spectra  of  Figures  6.6  to  6.8  have  been  sup- 
pressed to  provide  a  look  at  realizable  filters. 

The  probability  of  detection  for  the  threshold  detector 
was  shown  in  Part  IV  to  be 


Q,  =  SNR  M 


CO 

E 
k 


+  MQkBNR' 


Figures  6.17  through  6.20  show  Q,  as  a  function  of  SNR, 
BNR  and  M. 

The  probabilities  of  detection  for  point  targets  located 
at  specific  points  in  the  scan  have  also  been  plotted.   The 
probability  of  detection  for  a  specific  point  target  is  shown  as. the 


5o 


width  of  the  detector  is  varied.   As  width  is  increased,  so 
is  the  area  scanned  by  the  detector.   Figure  6.5  shows  the 
difference  in  spectra  between  a  rectangular  detector  and  a 
square  detector  which  have  approximately  the  same  surface 
area.   It  should  be  obvious  however  that  the  rectangular 
detector  scans  much  more  area  per  nutation  than  a  square 
detector  of  identical  surface  area. 

The  envelopes  of  signal  coefficients  (magnitude  only) 
for  a  point  target  are  plotted  in  Figures  6.6  through  6.8. 
Phase  information  depends  only  on  the  location  of  the  target 
relative  to  the  initial  point  for  the  nutation  cycle  which 
is  along  the  x  axis.   This  has  not  been  plotted.   Changes  in 
the  magnitude  due  to  changes  in  detector  size  are  shown. 

Figures  6.9  and  6.10  show  the  envelopes  of  coefficients 
for  the  signal  that  is  averaged  with  respect  to  a  Gaussian 
pointing  error.   The  bandwidth  of  the  averaged  signal  is 
much  less  than  a  point  target  as  expected.   The  figures  she 
the  spectra  for  two  different  standard  deviations  of  point- 
ing error. 

Matched  filters  (magnitude  only)  for  the  threshold 


Q'  =  SNR'i/M"  _- 


Z  f 
k 


( t  -\  f  *  

kK  oJ    k  i  +  MQkBNR^ 


I  |f  |  *  ± 

k   k    l+MQkBNR2 


where 


SNR'  -  h   \  an  ^7/1 

11  n=l       '  \J     2 
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Figures  6/21  through  6/24  show  Ql  as  a  function  of  SNR' , 
BNR  and  M. 

The  parameters  of  the  nutating  optical  system  have  been 
varied  to  show  their  effect  on  Q,  and  Q^. 
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VII.   CONCLUSIONS 

The  problem  investigated  in  this  thesis  was  to  determine 
an  optimum  processor  to  be  used  in  detecting  targets  with  an 
infrared  nutating  system.   The  nature  of  the  system  divided 
the  problem  into  two  areas:   optimizing  a  spatial  processor 
according  to  size  and  shape  and  optimizing  a  temporal  pro- 
cessor that  takes  the  output  time  varying  voltage  and  decides 
if  a  target  is  actually  present. 

General  equations  for  the  output  from  a  nutating  detector 
were  known  from  Samuelsson's  work,  however  these  equations 
had  not  been  solved  for  a  specific  detector.   At  the  same 
time  it  was  observed  that  the  form  of  the  covariance  function 
for  the  noise  offered  an  easy  solution  to  integral  equation 
of  the  Karhunen-Loeve  expansion  which  made  statistical 
detection  theory  appealing.   Some  work  had  been  done  on  th : 
form  of  temporal  processor  using  statistical  detection  the  iry 
but  this  was  limited  to  only  one  nutation  because  of  the 
common  background  noise  component  between  nutations. 

Harger's  derivations  provided  a  means  for  describing  a 
temporal  processor  using  statistical  detection  theory  that 
based  its  decisions  on  multiple  observations  and  included 
the  background  noise.   His  work  was  extended  here  to  include 
unknown  parameters,  amplitude  and  position,  which  more 
closely  characterized  the  system.   This  extension,  called  a 
threshold  detector,  was  shown  to  be  optimum  in  the  case  of 
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small  signal-to-noise  ratios  where  detection  is  most  diffi- 
cult.  The  derivation  also  used  the  Gaussian  assumption  to 
describe  the  background.   This  may  or  may  not  be  correct. 
But  in  practice,  the  detector's  performance  may  not  suffer 
greatly  if  this  assumption  is  wrong. 

To  implement  the  threshold  detector,  the  integrals 
describing  the  detecor  output  were  solved  using  the  Gaussian 
quadrature  method  of  numerical  integration.   Checking  the 
computed  coefficients  by  a  summation  provided  a  means  of 
accepting  the  validity  of  the  integrations.   The  form  of  the 
threshold  detection  system  using  a  rectangular  detector  was 
determined  and  the  frequency  spectrum  of  the  optimum  filter 
was  shown.   To  specify  performance,  the  probability  of 
detection  is  plotted  against  signal-to-noise  and  background- 
to-noise  ratios  and  the  number  of  nutations  on  which  the 
decision  was  based.   The  signal  spectrum  for  a  point  target 
and  background  power  spectral  density  is  also  plotted. 

A  designer  may  find  these  results  useful  in  developing  a 
system  or  investigating  the  performance  of  an  actual  system  as 
compared  to  the  optimum.   He  may  also  use  these  results  to 
determine  the  size  of  a  rectangular  spatial  filter  to  be  used 

Future  work  should  include  extending  these  results  to 
spatial  filters  that  are  circular  or  elliptic.   The  same 
numerical  integration  algorithm  could  probably  be  used.   The 
threshold  detector  was  found  to  be  highly  sensitive  to  back- 
ground correlation  lengths.   Thus,  some  form  of  adaptive 
processor  should  be  included  to  minimize  the  effect  of 
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background.   At  the  present  time,  a  measurement  program  is 
underway  at  Naval  Weapons  Center,  China  Lake,  California 
to  determine  average  background  correlation  lengths  but  for 
different  environments,  these  correlation  lengths  could  be 
expected  to  vary  considerably. 
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APPENDIX  A 
DERIVATION  OF  SIGNAL  SPECTRAL  COEFFICIENTS 

The  radiance  on  the  detector  is  periodic  because  of  the 
nutation  therefore  it  can  be  expressed  in  a  Fourier  series, 

°°     jnoo  t 
H(t)  =  E  Hn  e*   °  (1A) 

n 

where 

,  rT      -jnoo  t 
Hn=iJ   H(t)e     °   dt  (2A) 

o 

But, 

H(t)    =   /n' (r)T(r-p(t))    d2r    .  (3A) 

Substituting    (3A)    in    (2A) 

H      =    fd    r   N»(r)    ±1       dt    xfr-p(t)    e  °  (4A) 

n        J       ~  ~TJ0  ~~ 

The  transmittance ,  x,  may  be  expressed  in  terms  of  the 

inverse  Fourier  transform,  assuming  an  infinite  image  plane, 

c   9        -j27Tk-(r-p(t)) 
T(r-p(t))  =  JdZk  x*(k)e     ~     ~     .  (5A) 

Rewriting  the  exponent, 

k-(r-p(t))  =  k-r  -  pkxcoswot  -  pkysintoQt       (6A) 
Equation  (4A)  can  be  rewritten, 


H      =      d    r   N1 v 


T 


-j2   k-r 


j2ttAsiii(w    t+0)       -inco    t 
o  J      o 

e  e 


where 


■    xf^V? 


x  y 


tan"      kx/k    . 


Observe  that 

N'(k)  =  f, 


?  -  j  2  it  k  •  r 

d  r  N1 (r)e 


thus 


(7A) 


(8A) 


Hn  =  |d2k  N'(k)T*(l()ejne  £  /   dt 


T      -jn(io  t+9) 


jnAsinfw  t  +  6) 

j  v.  O 


C9A) 


Using  the  relation 


'n<*>    "   f/ 


-in<P      izsin$    ,. 
e    J         eJ  d$  , 


(10A) 


the  expression  for  the  coefficients  becomes 


H   =  fd2k  N'(k)x*(k)  e^nQ   J   2ttP  Jk2  +  k2 
n    J  ^-J       K~J  n I   K  \  x   y 


C11A) 


Let 


then 


<j)  =  tan    ky/kx 


-  n  -  j  n  d>    j  n  i 
j    e  J  Y  =  eJ 


Rewriting  (11A) 
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Hn    =    jn   Jd2k   N'(k)x*(k)    e-J^jlzvpJkl+k*     J.       (12A) 
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APPENDIX  B 
DERIVATION  OF  BACKGROUND  CORRELATION  COEFFICIENTS 


The  n   power  spectrum  component  of  the  background  at 
frequency  nw   is  given  by 


n  "  T  f  RB^ 


inco 


dT 


(IB) 


But  Rt,(t)  is  the  correlation  function  of  detector  output 
B 

averaged  with  respect  to  the  starting  time. 


Rb(t)  =  \  J      E[B(t  +  T)B(t)]dt 


(2B) 


Using  the  expression  for  B(t)  from  Part  II, 


Rb(t)  =  \j      dt  E 


/N,(r)i(r1-p(t+T))d2r1 


/N'(r2)T(r2-p(t)) 


d  ll 


(3B) 


Interchanging  the  order  of  expectation  and  integration, 


T 
Rb(t)  -  1/   dt  /d2ri  /d2r2  T^-fiCt+T)] 


,T[r2~£Ct)]E[N'c~i)N'(~2l] 

Assuming  stationary  background, 

E[N'(r1)N'(r2)]  =  ^(^-r^) 


(4B) 


(5B) 
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where  <j)'(r,-r~)  is  the  background  covariance  function  on  the 
image  plane. 

It  is  now  convenient  to  express  the  integrand  in  terms 
of  two-dimensional  spatial  Fourier  transform.   To  begin  with 


r  j27rk-r   ? 

x(r)  =  jT(k)e    ~  ~  dzk 


(6B) 


where  t  (k)  is  the  Fourier  transform  of  the  aperture  function 

r(r)  and  k  =  (k  , k  ) .   Then 

x   y 

T[r1-p(t  +  T)]r  [r2-p(t)J 

2  (  2  j27T(k    ■r1-k2T2) 

=   Jdzk1T(k1)    Jdzk2T*(k2)e  ~ L   -1   ~z    "L 

-j2Tr[k1    P(t  +  T)-k2p(-i)] 


(7B) 


Substituting  (5B)  and  (7B)  in  (4B)  and  interchanging  the 
order  of  integration, 

Rb(t)  =  /d2^  xCkj)  /d2k2  x*(k2) 

rl  -j2Tr[k1-p(t  +  T)-k2-p(T)] 

1  tJ   dt  e 


i  2 Trk  *  r 
•  /d2ri  e    -1  -1  /d2r2  ^(rrr2)e 

The  last  integral  becomes,  with  R  =  r-^  -  r2 , 


'*2vh'Z2 


(8B) 


Jd2R(J)'  (R)e 


-  j2TTk?-R- 


-  j  2Trk2  •  r. 


(9B) 


But  the  bracket  above  is  the  Fourier  transform  of   BM,  or 
the  Wiener  spectrum  of  background  on  the  image  plane.   Let 


-j2TTk2-R   2 
WB(k2)  =  J^B(?)e     ~Z    d  R 


(10B) 
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If  F_(k)  is  the  Fourier  transform  of  the  point  spread  func- 
tion and  WR(k)  is  the  Wiener  spectrum  of  background  on  the 
object  plane,  then 

wb(^  =  IV^I2   V^  (11B) 

This  relation  is  analogous  to  the  output  noise  spectrum 
expression  due  to  a  noise  input  to  a  linear  system.   In  terms 
of  W£(k)  the  last  line  of  (8B)  is 


W 
B 


(   j  2tt  (k-.  -k?)  ■r1   ? 
(k)  Je     -l    -L      -1  dZrx  =  W^(k)6(krk2)      (12B) 


Substituting    this    into    (8B)    and  performing    the    integration 
with   respect    to   k~ , 


RB 


(t)    =    /d2k    |x(k)|2    W£(k) 


,     fT      j2irk'  [p(t+x)-p(t)] 

-  J       e         "      ~  "  dt  (13B) 

o 


Since    p(t)    =    (p    cosco   t,  p  sinco   t)  , 


k    •     [p(t+i)  -P(t)  ]     =    kxpcoscoQ(t+T)  +  k    psinooo(t  +  T) 


-k    pcosoo    t    -    k    psinw    t 
x  o  yK  o 


where 


=    A(t)coso)    t    +    B(T)sinajQt    =    c  (x)  cos  (ojQt-ip  (t)) 


A(x)    =    p[kxcoso)oT    +   kysino)oi    -    kx] 
B(t)    =    p[-kxsinwoT    +    kycoscoQT    -    ky] 


C(t)    =     ^A2(t)+B2(t)    =    P  ^2(1    -    coscooT)(k2  +  k2) 
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and 


iKO  =  tan"1  *£l|  (14B) 


Therefore  the  last  integral  of  (13B)  is 

-.      rT       J2TTC    cos(co    T  +  lp) 

YJ   e  dt  =  J0(2ttc)  (15B) 


where  J  (')  is  the  zeroth  order  Bessel  function.   Therefore, 
the  correlation  function  RR(T)  is  given  by 

RB(T)  =  /lT(^|2  wb(^  Jo(47Tp^  sin  "T~)d2k     (16B) 

where 


k  =  Jk2  +  k2    and   d2k  -  dk  dk   . 
—    \  x   y  .,     ~x  ~y 

Finally,  the  power  spectrum  coefficient  8   is,  from  (IB) 


n 


6n  =  /|x(k)|2  W^(k)d2k 


,       rT  /  03T\-jnOJT 

i  J   J0  2a  sin  -§-  e     °  di         (17B) 


where  a  =  2"rrkp.   The  last  integral  becomes  with  x  =  w0T/2  , 


1  f^ 

—  I   J  (2a  sin  x)  cos2nx  dx 

T  J  o  ^         J 


1  (^ 
j  -  J   J  (2a  sin  x)sin  2n^  dx         (18B) 


2 
The  cosine  integral  equals  Jjj(a)  and  the  sine  integral 

vanishes  [2].   Therefore  on  using  (18B) 
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3n  =  JlT(^|2  I  FoC15^  I  2  WB(k)J^(27rpk)d2k  .       (19B) 

This  is  the  general  power  spectrum  expression  for  a  nutating 
detector  output.  The  power  at  frequency  noo  is  expressed  in 
terms  of  aperture  transform  i(k),  optical  transfer  function 

F  (k) ,  Wiener  spectrum  of  background  WR(k)  and  nutation 

2 
radius  through  J  (2iTpk)  .   The  integration  is  over  the  entire 

k-plane . 
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APPENDIX  C 
EXPECTATION  OF  A 


12 


The  essential  calculation  required  to  find 
A12((Zn))  =  EB[A12({Zn)|B)] 


(1C) 


for  Gaussian  background  B  is  the  expectation 


where 


I    =   E 


B 


exp|/      A(t)B(t)dt    -    c    /      B    (t)dt 
o  o 


c    =    M/N 


and 


2 


M 


1N   n=l 


C2C) 


B(t)  may  be  expanded  in  its  Karhunen-Loeve  (K-L)  repre 
sentation  provided  the  mean  square  value 


E(B  (t))  <  oo 


and  thus  the  integral  equation 


(3C) 


X.<j>.  (t)  =  /   RR(t,u)  *-  (u)du 
11      0 


0  <  t  <  T 


(4C) 


can  be  solved.   A.  are  called  the  eigenvalues  of  the  equation 

and  <t>    (t)  are  called  the  eigenfunctions . 
i 

Using  Mercer's  theorem,  the  correlation  function  of  the 
background  may  be  expressed  as 
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RB(t,u)  =  z    *k4>k(t)  <J.*k(u) 


(5C) 


The  K-L  expansion  of  B(t)  is 


B(t)  =  Z  bk<j)k(t) 


with 


T 
bk  =  /   B(t)  <J>£(t)dt 
o 


0  <  t  <  T 


(6C) 


where  the  coefficients,  b,  ,  are  unco rre la ted  and  because 
B(t)  is  Gaussian  the  coefficients  are  statistically  inde- 
pendent with  means  jj,  and  variances  A,  . 
Also  expanding  A(t) , 


A(t)  =  I  akc()k(t) 
k 


with 


ak  =  /   A(t)  <J>£(t)dt 


0  <  t  <  T 


(7C) 


Rewriting  in  terms  of  the  expansions 


I  =  E, 


exp 


oo  oo  oo  oo 

!      i\l    I   anbk  +  J1    I    anV  *kCt)**(t)  dt 
o     n  k  n  k 


c   /  Z    £  bnbf  <J>k(t)<|)JCt)  dt 
o   n  k 


(8C) 


Using 


6    =  /   4>  (t)d>*(t)dt 
xy       yx^  J^y^   J 
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I  =  E-. 


bn 


exp 


I    (i  a  b*  +  \   a*b  )  -  c  I    b  b* 
2   n  n    2   n  rr       n  n 
n  n 


(9C) 


Using  the  independence  of  the  coefficients 


I    =    IT    E,       exp 


n 


bn 


x  a   b*    +    75-  a*b      -    c 
2      n   n         2      n   n 


n 


CO 

=   n   exp 
n 


n 


4c 


Ei       exp 
bn        r 


-c 


a      2 
n 


n         2c 


(IOC) 


For  Gaussian  random  variables 


E(wxx*)  = 


exp 


w  m  m*/ (l-2wA) 
x  x  *•      J 


1  -  2wA 


where 


m  =    E(x) 
x    v  J 


X      =  E[ (x-m  )  (x*-m*)] 


(11C) 


and  w  is  a  constant 


Let 


x  =  b   -  a  /2c 
n    n' 


and 


m   =  u   -  a  /2c 
x    n    n 


then 


E,   exp 
bn   r 


1  a  2 

u  n 

-CD  "  o— 

n  2c 


=  exp{-clun-an/2c!2(l  +  2cXij} 


(12c) 


(1  +  2cX  ) 
^       n  ■ 


Substituting  (12C)  into  (IOC)  ,  and  separating  terms 
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B  .  I      »      |a    |2X 

I   -B    CI   ♦    2CXJ"1   expUl   r^W1 
n  n  n 


exp 


od     y  *  (  y     -a     /  C  ) 

2  1    +    2cA 

n  n 


exp 


00  y    fy*-a*/c") 

c_  y  MnVMn      ir     ; 

2    „  1    +    2cX 

n  n 


(13C) 


This    result   can   be   written    in   closed    form  by   using    the 
definition  A      =    fa,  <b    ) 


oo      a         a  T  T 

1  Z    1  +n2cA      =    m  f      dtlA^tl)    /      dt2A(t2)q(t      t    )       (14C) 
n  n  o  o 


where 


q(t1,t2) 


2M 

N 


1  +  2cX      Yn v    r Ynv    2J 
n  n 


(15C) 


is    such    that    operating    on    q(t,,t~)    with 


T    T 
/    / 
o    o 


M~ 


(•)dt1dt2 


(16C) 


and  using  Mercer's  theorem  yields 


N 


^q(t1,t2)  +  /   RB(T2,t3)  q(t3,t1)dt3  =  RBCT1,t2)   (17C) 


Likewise , 


V*  (\i    -a  /c)  +  y  fy*-a*/c) 
n  *•  n   n       n   n   n 


n 


1  +  2cX 


n 


\   /  (mb(t)  -  ^-)h(t)  dt 


C18C) 
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where 


2y 


n 


h(t)  =  E  ..  M     "    (f)  (t) 
^  J  1  +  2cX  ^nK    J 

n       n 


(19C) 


such  that 


Nh(t)  +  2M  /   RB(t3,t)  h(t3)dt3  =  MmbCt) 
o 


(20C) 


Taking  the  logarithm  of  both  sides  of  equation  (13C)  and 
using  the  equations  above 


in   I  =  Z  (1  +  2cAn) 
n 


-1 


N 


T 


+  -rvr  /    dt, 
4M  1 

o 


I  j^W  -  W) 


T 
/   dt 
o 


2   M 

£  E   (Z  (t~)  -  S  (t9)) 

m=l 


qCt1,t2] 


/ 

o 


mb(t) 


N 
M 


M 


g.  z  (z£(t)  -  szM) 

X/   -L 


h(t)dt 


In   I    =  Z  (1  +  2bXR) 
n 


,   T  T 


M         M  M  M 

£=1       m=l         £=1        m=l 


M 


M 


M 


m 


E  S  ft  )   E  Z  (t  )  +   E  S  (t  )  E  S  (t,) 
=  1  m  L      1=1    *   -1    £=1     i  m=l  m  ^ 


q(t1,t2)dt1dt2 


/   m  (t)h(t)dt 
o 

T  M 
2  Z1   E  (Z  (t) 
M  o  1=1 


S£(t))  h(t)dt 


(21C) 
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The  likelihood  ratio  desired  is 


A12({Zn})  =  exp 


7      M   T 

£  Z   /   Z  (t)S  (t)dt 
N  n=1  0 


1   M   T   2 
n=l  o 


•  I 


(22C) 


Taking  the  logarithm  of  both  sides  and  substituting  the 
expression  for  I    (I), 


n 


M        T 


£n   Al2({Zn})    =   fi     \    f      Zn(t)Sn(t)dt 

n=l    o 


1      M        T      2 

--If  S^(t)dt 
N         i  n^    ^ 

IN   n=l    o 


+  I  "n  C1  +  f  V1 


T  T    ,    M 


+  m'  ! 

O    O 


M 


£=1  m=l 


MM  MM 

1  =  1  m=l  m=l  1=1 


M  M  . 

\    W     'S      Sm(t2)       q(t1,t2)dt1dt2 
,=  1  m=l 


T  ?      T     M 

/      mb(t)h(t)dt    +    §  f        Z       (Z£(t)  -  S£(t))h(t)dt 

O  O        x.—  1 

(23C) 
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APPENDIX  D 

KARHUNEN-LOEVE  REPRESENTATION  OF  EQUIVALENT 
SIGNAL-TO-NOISE  RATIO 


The  variance  of  the  threshold  statistic 


M   rT 


W  =   I   /   Z  (t)q(t)dt 
n=l  o 


where 


q(t)  -  ± 


f(t)  -  /   £(u 
o 


)q(u,t)du 


(ID) 


under  hypotheses  H   and  H,  is 

o      1 


VAR 


T 

=  ^/   q2(t)dt  +  M2  fdfq(t)  /duq(u)Rft,u)   (2D) 


Substituting  (ID)  into  the  first  term  of  the  variance 
expression 


NM  f      -  2  2M  f  -?  "     fT  - 

^J   qZ(t)dt  =  ^J   dt  fZ(t)  -  2f(t)  J   f(u)q(t,u)d 


u 


+J  J       f(u)f (v)q(t,u)q(t,v)dudv 
o  o 


(3D) 


Using  the  Karhunen-Leeve  representation 


NM 
z 


/v 


r^A  +  2M 

(t)dt    =    -^- 


r 


k  I     K   l     o       K         l 


2    n    f,  f ; 


2M    . 

T  A 


k    £    ,    t 2M 

k      £  1+-TT     1  o 

N      n 


T  T 

—  /      ^k(t)^*(t)dt  /    d>    (u)d>*(iOdu 
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-  oo  co  co  _  .    2m  A      2m  A      T 

+  Z  Z  Z  Z  f  f*  N   9mm N   n       ({,  (t)4)*(t)dt 

v  0       k  £  1   2m  .   ,  J  2m  ,   J    mv  ;  n^  ; 
k  I    m  n        1  +  -rr-  A   1  +  t=-  a   o 

N   m     N   n 


T  /-T 

/   4>k(uH*(u)dt  /   <f>nO)4>*0)d 


(4D) 


Using 


6nk  "/   ♦„Ct)*{|(t)dt, 


then 


NM 


■T 


/   42(t) 


di    2M  "  ~    ~      I  N   k    + 

\      .+  2M  ^ 

.2 

N  t   kk  I ,   2M 

N   k, 


2M  , 
IT  Ak 
1+  2M 

1  TT  Ak 


(5D) 


The  second  term  may  be  expressed  in  a  similar  manner  to 


give 


—  E    Vk  ( 2T1  — 

N     k  \l  +  f  Xk 


*  (™  v  i 


i 


VN;      k    *kk II7ZH   A 


Ak 


N       k/ 


VAR  =  f  J   VJ  — 4m— 

k  1+ir  Ak 


co) 


The  equivalent  signal  to  noise  ratio  in  Part  II  is 


M 
Z 

n=l 


-T  . 


in        r  i 

Z  ■  an  J   £(t)q(t)dt 


(7D) 


^  VAR 


The  K-L  representation  of  the  numerator  is 


M  rT    _  M 

Z      a  f(t)q(t)dt   =      Z   a 

n=l  o  n_J- 


I  J    lfkl2/(1  +  f  ^k 


.     (8D) 
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d  may  now  be  represented  by  dividing  (8D)  by  the  square  root 
of  (6D) 


d   = 


0  M  rli     |2/m     L  2M    ,     > 

1  Z      an   ^lfkl    ^1+~W  Xk) 
n=l 


2M       .-     ,2,        +2M 
N    M±k|    /    1      N    Akj 


d  = 


M 

«      E      a      M 
M        ,      n 
n=l 


E|fk|2/d+f   Ak) 


v7!^" 


(9D) 
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APPENDIX  E 
NON-ZERO  TERMS  FOR  NUMERICAL  INTEGRATION 

Equation  (5.5)  for  the  radiance  function  for  a  point 
target  may  be  written  as 

.2     ,        2       2       -j (p    +Q    )       .  .     . 

T_2_     f-x-y  Jvxxy^    s  max    smgy 

,2    Je  e  x  y 

•e"^   jfc^x2  +  y2     )dxdy 


where 


and 


a  =  w//2~a 
6  =  h//2a 
c   =    /2~p/a 


^2    r  s 

P  =   —    (r    -x   ) 
a         x      o 


Q     =    {1     (r     -y     ) 

x         a         y      o 


<(>   =    tan"1   y/x  (IE) 


Let  E(x,y)  be  that  part  of  the  integral  that  is  even  about 
both  planar  axis  which  is 

Using  Euler's  relation  the  integral  becomes 
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I  =  JE(x,y)  [cos  (Px+Qy)  +  j  sin  (Px+Qy)  ]  [cosn<J>  -  jsin  n<J)]dxdy    (3E) 
Using  trigonometric  identities 

r 

jE(x,y)  _Z   T.  dxdy 


I  = 

i=l 


where 


T,  =  cos    Px    cos    Qy    cos    n<J> 

T?  =  -sinPx   sin  Qy    cos   n<$> 

T_  =  -sin   Px   cos   Qy   sin  ncf> 

T.  =  -cos    Px   sin  Qy   sin  n<j> 

Tr  =  -j    sin   Px   cos    Qy   cos    ncf> 

TV  =  -j    cos    Px   sin  Qy   cos    n<$> 

T-  =  -j    cos    Px   cos    Qy   sin   n<J> 

T8  =  +  j    sin   Px   sin  Qy    sin   n<f>                                                     (4E) 

To  save  computation  time  it  is  necessary  to  determine 
which  of  the  eight  integrals  are  always  zero  and  for  which  n 
they  are  zero. 

Observe  that 


cos   n<fr   -   R     ejn<l>   =    R    (e^)" 


e  e 

n 

=    R    (coscf>      +    j    sin   40 


x  .         y 

+  J  — — 


n 


2.2  /    2.     2 


x    +y  yjx   +y 

Nn/2 

R„    (x  +  jy)n.  (5E) 


2M    2  e 

.x   +y    / 
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Likewise 


n/2 


sin  nty 


2A  2 
,x  +y 


I  (x+jy) 

m  V.     J  J   J 


n 


(6E) 


The  first  term  in  brackets  is  always  even  and  may  be  made 
part  of  E(x,y).   The  second  term  in  brackets  may  be  expanded 
by  the  binomial  expansion  with  the  general  term 


,  * n-k  ,  .  . k 
(x)    (jy)   . 


The  first  integral  of  (4E)  becomes 


(7E) 


I,  =  IE(x,y)cos  Px  cos  Qy  R  (x+jy)   dxdy 


(8E) 


n 


For  I,  to  be  non-zero,  Re(x+jy)   must  be  even  in  x  and  even 
in  y.   From  (7E)  this  is  true  only  of  even  n;  for  odd  n,   I, 
will  always  be  zero. 

Examining  the  other  seven  integrals  yields 

jE(x,y)[cos  Px  cos  Qy  cos  n0 

+    j    sin   Px   sin   Qy   sin  n<f>]dxdy        n    even 

jE(x,y)[cos    Px   sin   Qy   sin  n<J> 

+    j    sin   Px   cos   Qy   cos   n<£]dxdy        n   odd 

(9E) 
Observing  that  symmetry  about  the  x-  and  y-axis  is  necessary 
for  this  integral  to  be  non-zero  and  integrating  only  when 
the  integral  is  non-zero  permits  one  to  integrate  only  over 
the  first  quadrant. 


107 


APPENDIX  F 
SUMMATION  CHECK 


The  signal  for  a  point  target  was  shown  in  Part  II  to  be 

H(o)ot)  =  /N'(r)T(r-p(t))  d2r  .  (IF) 

Assuming  an  infinite  image  plane,  H(w  t)  may  be  written  in 
terms  of  the  Fourier  transform  as 

f  -j2iTk-p(t)   ? 

H(o)  t)  =  J  N'  (k)x*(k)e     ~  ~     dzk  (2F) 

For   a    rectangular   detector    (2F)    may   be   written    as 

/-2tt   a    fk    +k    )       .         ,  .      ,, 

K   x      yJ    smTrwkx  smirhky 

TTK  7TK 

x  y 


f  j2irk-(p(t)-r  )  -j27r(kxxo+k   y   )    d2k      (5F) 

le 


•    e 


Let 

x   =    /2    airkx. 
y   =    /7  airky 
and  noting    the    integral   may   be    separated   into    twoparts 

/2 

,     ,        2  i —  x(pcoso)    t  -   (r    -x    )) 

H(u>    t)    =   I    fe"x      sin  ™L_  e'   a  °  x      °        dx 

0  ^  j  /2a 

hy 


sm 


..     r        2             /2a      j —  yfpsinoj    t  -  (r    -y    )) 
I      e-y e     a  °  y      °      dy 

(4F) 
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Since  the  integration  is  -°°  to  +°°  only  even  integrands  will 
be  non-zero.   Rewriting  (4F) 


H(co    t) 

^    o    J 


T 


/ 


sm 


wx 


/2a  /2       . 

cos    —  x[pcosco    t 

x       \o  o 


(rx-xo)]  dx 


2-f 


hy 

in  — L- 


e"y — —  cosf—  y[psinwot  -  (ry-yQ)] 


The  general  form  of  the  integral  is 


(5F) 


2 


R  =  R 


CBD  -  -/  e"U  ^ 


-^  cos  (cu)  du 


(6F) 


Differentiating  with  respect  to  B 


dR(B)    2  f         -uz     „  .   , 

,p  ^  =  —  J   e     cos  By  cos  Cu  du 


(7F) 


But  the  solution  to  (7F)  is  known  [2]: 


dR(B) 
dB 


2/fT 


(B-CT       (B+C) 


+  e 


(8F) 


Since  the  constant  of  integration  R(0)  =  0, 


R(B)  = 


2/? 


/ 


(v-c) 


f 


_    (w+c) 


dv  +  J 
o 


dw 


(9F) 


Let 


P  = 


v-c 


and 


Q  = 


w+c 
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then 


R(B)  = 


I, 


B-C 
A-f  2   e"p2  dP  + 


riTJ-C/2 


B+C 
—  f   2 


Erf  (^)  +  Erf  (^) 


2     n2 
e"Q   dQ 


(10F) 


where 


Erf(z)  =  —  f 


rz  2 

e  }      dy  . 


(11F) 


Substituting  the  result  (10F)  into  (5F) 


HO  t)  =  4- 

v  o  J         4 


Erf 


T 


pcosw  t  -  r  +x 
o     x   o 


7a 


w 


"4 


pcosa)  t  -  r  +x 

O       X   o  J 


Ta 


w 


+  Erf 


-Erf 


/2 


ps  in 


w  t  -  r  +y 
o     y 


•l-i') 


psinoj  t  -  r  +y 
o     y  7  o 


7a 


(12}  j 


The  mean  square  voltage  of  the  background  may  also  be 
solved  in  a  like  manner. 

From  Part  V,  the  mean  square  voltage  was  shown  to  be 


-2 
eB 


||T(k) 


2  |FQ(k)|2  Wfi(k)  d2k 


(13F) 


which    for   a    rectangular    detector   becomes 

2a2,  2 


.sin2  (irk   W) 
2  0      /  K      x    J       e 

eB  =  V2a/-      TX 


■/ 


4TT      "k' 


77     k' 


a     +  (2Trk    ) 
x 


T  dk 
2  x 
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•aB-2B 


2 
sin  (irk  h) 
y  J 


/sin 


V 


a    2  2,  2 

-4tt  a  k 

e       y 
— ~ *  dk 

ez+(2Trky)z     y 


Letting 


t  =  2-nk 


the  general  form  of  the  integral  to  be  solved  is 


(14F) 


R 


oo  -  i  / 


9  2  +  2 

00   ■  2 ,  .  /0.>    -at 

sin  (y t / 2 j   e 


t2/4 


2  .  -2 
a  +  t 


C15F) 


Taking  the  first  and  second  derivatives  one  has 


dR(y) 
du 


2-2 
a  t 


1  f   sin(2ut)   e , 

o      L      a  +t 


and 


(16F) 


d2R 


R(y)     2  f°°  cos(2yt)   -a2t2 

^2      tt|    2i+2    e 
du        0   a  +  t 


(17F) 


The  solution  to  the  second  derivative  is  known  [2] 


~         2  2 

d  R(y)  _  e 

,  2       2a 
du 


-ay  r    r      r  V    -\ 

e    Erfc  (aa-  -S-- J 

^     2a^ 


+  eay  Erfc  (aa  +  ^-)  ] 


where 


(18F) 


/ 


y 


Erfc(z)  =  —  J   e  '   dy 


(19F) 


is  the  complimentary  error  function. 

With  lengthy  and  tedious  calculations  it  can  be  shown 
by  using  the  relation 
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R 


(y)  =  J   dv  J   F(w) 


dw 


(20F) 


with  initial  conditions  equal  to  zero  and 


F(w)  -  dRM 
du 


(21F) 


that 


a 


1  -  Erfc  (H-) 


2a 


1  -  e 


-(y/2a) 


/ttu 


2  2 
a  a 


2a' 


ay  r  r   /     u  % 
e    Erfc  (aa  +  •=— ) 

v      2a' 


-ay  ~  r    ■  r  y  n 

+  e   p  Erfc  (aa  -  4y — 1 


2  Erfc  (aa) 


The  mean  square  voltage  thus  becomes 


(22F) 


el   =  Aafrol    R(w;a)R(h;6)  . 


B 


B 


(231 
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APPENDIX   H 
C 
C 

C  COMPJTER    PROGRAM    TO    CALCULATE    THE     FOURIER    COEFFICIENTS 

C  OF    THE     BACKGROUND    CORRELATION    FUNCTION, A     POINT    TARGET 

C  AND    AM     AVERAGED    TARGET    OVER    A    GAUSSIAN    POINTING    ERROR 

C  FOR    A    NUTATING    OPTICAL    SYSTEM    WITH    A    RECTANGULAR 

C  DETECTOR. 

C 
C 

c 

C  AA    =    BACKGROUND    CORRELATION    FUNCTION    COEFF(ZERO    ORDER) 

C  ALFA    =     BACKGROUND    CORRELATION    COEFFICIENTS     (VECTOR) 

C  FNO    =    SIGNAL     FOURIER    COEFFICIENT     (ZERO    ORDER) 

C  FN    =    SIGNAL    FOURIER    COEFFICIENTS     (VECTOR) 

C  FAVGO     =    AVERAGED    SIGNAL    COEFF    (ZERO    ORDER) 

C  FAVG    =    AVERAGED    SIGNAL    FOURIER    COEFF    (VECTOR) 

C  RHO    =     MUTATION!     RADIUS 

C  SIG    =    STANDARD    DEVIATION    OF    BLUR    CIRCLE 

C  WT     =    WIDTH    OF     DETECTOR     I X-DI RECT I  ON) 

C  H    =    HEIGHT    OF    DETECTOR     ( Y-D IP  EC T  ION ) 

C  AL    =    INVERSE    CORR    LENGTH    OF    BACKGROUND    ( X- D IR ECT I  ON ) 

C  BE    =     INVERSE    CD^R    LENGTH    OF    bACKGROUND     ( Y-DI RECT I  ON) 

C  PO    =    STANDARD    DEVISTION    OF    POINTING    ERROR 

C  XO    =    DETECTOR    COORDINATE    X-D I RECT I  ON 

C  YO    =     DETECTOR    COORDINATE     Y-OIRECTION 

C  RX    =    POINT    TARGET     COORDINATE    X-DIRECTION 

C  RY    =    POINT    TARGET    COORDINATE    Y-DIRECTION 

C  R.RANGLE    p    POLAR    COORDINATES    OF     POINT    TARGET 

C  NS    =    NUMBER    OF    FIRST     COEFFICIENT     DESIRED 

C  NUM    =    NUMBER    OF     LAST    COEFFICIENT    DESIRED 

C  IAVE=1       "SIGNAL"     WILL    CALCULATE    AVERAGE     SIGNAL    COEFF 

C  IAVE=0       "SIGNAL"     WILL    CALCULATE     POINT    TARGET    COEFF 

C 

C 

IMPLICIT  REALMS  (A-H.O-Z) 

DIMFNSICN  ALFA(120) 

C0MPLEX-M6    FN(  120)  ,  FNO,  FAVG(  120)  , FAVGO 

COMMON    /DTPRM/     RHO , S I G, W I , H, AL , BE , RX , RY , XO , YO 

COMMON     /AVPRM/     MUM, I  AVE, PO 

COMMON     /START/     NS 

COMMON    /ANGL/     RANGL 

RFAD    2, RHO ,SIG, WI ,H,AL ,BE 

READ    2  t XO. YO, R, ANGL 

READ    3,NS,NUM,I AVE ,P0 

RANGL  =  0  .01  745DOANGL 

RX=R*DCGS ( RANGL) 

RY  =  R>"DSIN(  RANGL) 

RANGL=DATAN?(RY-YO,RX-XO) 

CALL    SUMCK 

CALL    BKGRN(ALFA,AA) 

CALL     SIGNAL! FN. FNO) 

WRITE     (9)     RH0,SIG,WI,H, AL, BE.RX,  RY,XO,  YO,  PO 

IAVE=1 

IF( IAV E.EO. 1)     RX=O.ODO 

IF( IAVE.EO.l )     RY=O.ODO 

CALL    SUMCK 

CALL    S  IGNAL( FAVG, FAVGO) 

WRITE     (9)     FNO, FN, AA, ALFA 

WRITE  I  9)  FAVGO, FAVG 

1  FORMAT (10  15) 

2  FORMA T( 8D1 0.0) 

3  FORMAT ( 31 5,  ID  10.0) 

4  FORMAT (8E10. 4) 
STOP 

END 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 
c 


c 

c 
c 
c 


SUBROUTINE  BKGRN(VP,V) 


SUBROUTINE  BKGRN  COMPUTES  THE  COEFFICIENTS  OF  THE 
BACKGROUND  CORRELATION  FUNCTION  UP  TO  120  ORDERS. 
HIGHER  ORDERS  MAY  EASILY  3E  COMPUTED  BY  ONLY  CHANGING 
THE  DIMENSION  STATEMENTS. 

ZP  =  VALUES  OF  RADII  FOR  DIFFERENT  ANNUL  I 

Z  =  WORKING  VECTOR  FOR  RADII 

MAXR  =  MAXIMUM  RADIUS  OF  OUTER  ANNULUS 

NIZ1  =  DIMENSION  OF  ZP 

VP  =  WORKING  VECTOR  IN  FIRST  PART  AND  RETURN  VECTOR 

IN  THE  SECOND  PART 

V  =  ZERO  ORDER  COEFFICIENT 

IMPLICIT    REAL-'-S     (A-H,0-Z) 

DIMENSION    VP< 1) ,DN(50, 120) 

DIMENSION    DI( 50),ZP( 20) ,Z(50) 

DIMENSION    R(24) , W(24) ,TH(25 ) 

COMMON    /AVPRM/    NUM,IAVE,PO 

COMMON    /DTPRM/    RHO , SI G, W I ,  H, AL  ,  BE , RX, R Y , XO, YO 

COMMON    /BKPRM/    A,B,C,D 

COMMON     /NFIRST/    NFIRST 

DATA    MAXR/4/,NR/l/,NIZl/20/ 

DATA    Z°/0.  00  0,1 .200,2. 4 DO, 5. 5 DO ,8. 6 DO ,11 .800,15.00,20 
ID 0,25. 0  00,3  0. 00 0,40.00  0,50.03  0, 60.0C0, 70. 0D0,80.0D0, 
290.0D0, 100.000, 150.0  00,20  0.0  00, 40  0.0  00/ 

NFIRST=0 

PI    =    DARCOS(-l.DO) 


CALCULATE  CONSTANTS  FOR  INTEGRAL 

CONST    =    AL*BE*SIG**2/PI/PI 

A=WI/2.D0/S  IG 

B=H/2.D0/SIG 

C  =  (AL*SIG)**2 

D  =  (BE*SIG)**2 

E  =  RHO/SIG 


C 
C 
C 
C 

c 
c 


SET  ANNULI  FOP  SPECIFIC  PARAMETERS 

DO    150     I=1,NIZ1 
150       Z(I)=ZP(I):/E 

I    =    0 

IFCZtNI Zl) .GT  ,MAXR)GO    TO    210 
DZ    =    PI /NR/E 
MMM=50-NIZ1 
DO    200     1  =  1,  MM 
ZU+NIZ1)     =    Z(NIZ1)    +    n-oz 
IF(Z( I+NIZ1) .GT.MAXRJGO    TO    210 
00    CONTINUE 


C 

c 

C 

C 


MAX    =    NUMBER    OF    ANNULI 
10       MAX=I+NIZ1-1 

PRINT    4,RH0,SIG,H, WI , AL,BE, A,B,C,D,E,DZ 
PRINT    4, ( Z( I  )  , 1=  1,MAX) 
DO    300    1=1, MAX 

COMPUTE    RADII , WEIGHTS     AND    ANGLES    FOR    EACH    ANNULUS 
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c 
c 


c 
c 
c 
c 
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120 


C 
C 
C 

c 


270 


C 

C 

c 
c 


130 
280 
300 


400 


600 


C 
C 
C 


500 


700 


1 

2 

3 

4 

12 

13 

14 

16 


CALL  WAA(Zd)  ,Z(  1  +  1)  ,R,TH,W,NPQ) 


NFIRST=  1 
Did)  =  O.DO 
DO  115  NN=1,NUM 
DN(I,NN)  =  O.DO 


COMPUTE  FUNCTION  FOR  EACH  RADIUS 

DO    280    K=l,NPQ 

CALL    3ESSL     ( R { K ) * E , VP , VJN , NUM) 

RV  =  VJN-  v2/DEXP(P (K  J**2) 

EP=OEXP(R( K)**2) 

DO    12  0    NN=1,NUM 

VP(NN)=VP(  NN)**2/EP 

REC=0.D0 


COMPUTE    FUNCTION    FOR    EACH    ANGLE 

DO  270  J=1,NPQ 

REC=REC+FU(R(K)  tDCOS(THU)  )  ,  R(  K  )  *- DS  IN  (T  H  (  J  )  )  ) 

V=O.DO 

DO  130   NN=ltNUM 

DN(  I  ,NN)  =DN(  I  ,NN)  +  REC<:VP(  NN)*W(  K) 

DI(I)  =  Did)  +  REC*RV*W(K) 

CONTINUE 


SUM  VALUE  FROM  EACH  ANNULUS 


DO 
J  = 

V  = 

v=v 

PRI 

NBQ 

PRI 

PRI 

DO 

VQM 

VP( 

DO 

J  =  M 

VP( 

IF( 

CON 

VP( 

PRI 

PRI 

CON 


400  1=1, MAX 

MAX  +1-1 

V  +  D I  (  J  ) 
*CONST 

NT  3  ,N,V,D1(MAX  ) 
=  0 

NT  13,NBQ,V 
NT  4, (DI( I )t 1=1 ,  MAX) 
500  NP=1,NUM 
=  0.DO 
NP) =0. DO 
60  0  1=1,  MAX 
AX-d-I 

NP)=VP(NP)+DN( J,NP) 
DN( J,NP) .  GT. VQM)  VQM=DN( J,NP) 
TINUE 

NP)=VP(NP)*CONST 
NT  13,NP,VP(MP) 
NT  14, NP,VQM, DN(MAX,NP) 
TINUE 


CALCULATE  SUM  FOR  INTEGRATION  CHECK 

SUM  =  V 

DO  700  J=1,NUM 

I=NUM+1-J 

SUM=SUM  +  2.D0WP(  I  ) 

CONTINUE 

PRINT  12, SUM 

RETURN 

FORMAT ( 1015) 

FCRMAT(  BE  10.0) 

FORMAT ( 2110, 1P2E15.6) 

FCRMAT( 1F8E15.6) 

FORMAT ( 5X,G2  0.13) 

FORMAT 15X,  14, 6X, G20  .13  ) 

FORMA T( 5X,I4,6X,G20.13,6X,G20.13) 

FORMAT (3( 1PD20. 12)  ) 

END 
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c 
c 
c 
c 
c 
c 
c 


FUNCTION    FU(X,Y) 


FUNCTION    FU     IS     AN    AUXILIARY    R3UTINE    FOR     BKGRN    THAT 
COMPUTES    THAT    PART    OF    THE    INTEGRAND    NOT     RADIALLY 
SYMMETRIC. 


IMPLICIT    REALMS     (A-H,0-Z) 

COMMON    /BKPRM/    A,B,CtO 

Pl=l.DO 

IF(X.NE.O.DO)     P1=DSIN(A*X)/(A*X) 

IF(Y.NE.O.DO)     Pl=Pl*DSIN(B*Y)/(  B^  Y  ) 

FU=P1**2/(C+X**2)  /(D+Y**21 

RETURN 

END 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 
c 


c 

c 
c 
c 


150 


SUBROUTINE  SI GNA L( REC , RECO) 


SUBROUTINE  SIGNAL  COMPUTES  THE  FOURIER  COEFFICIENTS  OF 
A  POINT  TARGET  3R  AVERAGED  COEFFICIENTS  OVER  A 
GAUSSIAN  POINTING  ERROR 

ZP  =  VALUES  OF  RADII  FOR  DIFFERENT  ANNULI 

Z  =  WORKING  VECTOR  FOR  RADII 

MAXR  =  MAXIMUM  RADIUS  GF  OUTER  ANNULUS 

NIZ1  =  DIMENSION  OF  ZP 

REC  =  RETURN  VECTOR  OF  COEFFICIENTS 

RECO  =     ZERO  ORDER  COEFFICIENT 

IMPLICIT  RE4L^8  (A-H,0-Z) 

DIMENSION  3N(120) 

DIMENSION    DI (50) , DNJ( 50,120 ) , DN ( 50, 120 ) , Z ( 50 ) , ZP ( 20 ) 


i  <;  t  r 


ON    R<24) ,  W( 24) ,TH{ 24) 


RECO, VS,F,WP 


COMMON 
COMMON 
COMMON 
COMMON 
COMMON 


RX,RY,XO, YO 


5C0 ,S.6D0,1 1 .800,15 .00,20 


C0MPLEX'U6    PEC(1  ),VEC(120), 
COMMON     /£VPRM/    NUM,IAVE,PO 

/DTPPM/    RHC,SIG,WI,H,AL, BE 

/SIPRM/    A,B,C,PX,QX 

/NFIRST/    NFIRST 

/START/    NS 

/ANGL/     RANGL 
DATA    ZP/O.ODOtl- 200,2.400, 
100,25.000, 30.000, 40 .000, 50. 000, 60. 000, 70. 000, 80. 00 0, 
290.00  0, 10  0.000,150.000,200.000,400.000/ 
DATA    MAXR/7/,NR/l/,NIZl/20/ 
NFIRST=0 

VS  =  DCMPLX( 0.00, -1.00) 
PI    =    DAPCOS(-l.DO) 
PI2    =    PI/2. CO 
KKK=0 


CALCULATE    CONSTANTS    FOR    INTEGRAL 


2  +  PO** 2) 


200 


GAMA=SIG 

IF(  IAVE .EC.l)     GAMA  =  DSORT(SI  G- 

SG   =    DSQRT(2.D0KGAMA 

A=WI/SG 

B=H/SG 

CCNST=A*B/P1/PI 

C    =    l.DO 

PX=-XO 

QX=-YO 

IF(IAVE.EQ.O)     PX=PX+RX 

IF(IAVE.EO.O)     3X=QX+RY 

PX=PX?DSQRT(2 .DO)/GA^A 

QX=QX*DSQRT(2.D0) /GAMA 

E    =    2.D0-PH0/SG 


SET    ANNULI     FOR    SPECIFIC    PARAMETERS 

DZ    =    0.00 

DO    150    I=1,NIZ1 

Z(I)  =  ZP(I  )/(2.D0:'E) 

I     =    0 

IFU(NIZl)  .GT.MAXRJGO    TO    210 

DZ=15.D0*PI/E 

MMM  =  50-M  Zl 

00    200    I=1,MMM 

Zd+NIZl)     =    Z(NIZl)     +    I*DZ 

I F(Z( I+NI Zl) . GT.MAXRJGO    TO    210 

CONTINUE 

PRINT    17,ZII+NIZ1) 
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c 
c 

C      MAX  =  NUMBER  OF  ANNUL  I 

210  MAX=I+NIZ1-1 
GO  TO  212 

211  MAX=I 

212  PRINT  4  ,RHC,SIG,WI  ,H,AL,BE, RX,RY,XO,YO, PX,QX 
215   PRINT  4, (Z( I ) , 1=1, MAX) 

220   DO  300  1=1, MAX 
C 

c 

C  COMPUTE    RADII , WEIGHTS     AND    ANGLES    FOR    EACH    ANNULUS 

C 

CALL    W£A(Z(!),Z(I+1),R,TH,W,NPQ) 

NFIRST=1 

NPR=NPQ 

DI(I)    =    O.DO 

DO    115     NN=NS, NUM 

DNJ( I,NN)=O.DO 
115       DN(  I,NN  )    =    O.DO 

IHKKK.  EO.l  )     GO    TO    300 

DO    280    K=1,NPQ 

RK2=R(K)**2 

IF(RK2.GT.174)     GO    TO    285 

EP=DEXP (RK2) 

IF(EP.GT.1.D    50)     GO    TO    285 

CALL    BESSL     ( R ( K ) * E , BN T VJN , NUM) 

RV=VJN/ EP 

PECO=F( R(K) ,TH,0) 
C 

c 

C      COMPUTE  FUNCTION  FOR  EACH  RADIUS 
C 

RECO=RECO*RV*W(K) 

PP=RECO 

DI( I )=DI( I )+PP 

DO  130  NN=NS,NUV 

BN(NN) =BN(NN)/EP 

REX(NN) =F(R(K) , TH , NN) 

REC(NN)=REC(NN)*BN(NN)*W(K) 

PP=REC( NN) 

DN( I,NN)=DN( I ,NN)+PP 

PP=REC(NN)  •:  VS 

DNJ( I ,  NN)  =  DNJU,NN)+PP 
130   CONTINUE 
280   CONTINUE 
300   CONTINUE 
C 
C 
C      SUM  VALUE  FROM  EACH  ANNULUS 


C 


VZ    =    O.DO 
VI    =    O.DO 
V2    =    O.DO 
DO    400     1=1, MAX 
J    =    MAX    +1-1 

IF(DABS(DI(J)).3T.CABS(VZ))     VZ=DI(J) 
I  F  (  D I  (  J  )  .LE.O.DO  )    VI    =     Vl  +  DKJ) 
400       IF  (DI  (  J)  .GT.O.DO)     V2    =    V2+DKJ) 
RECO=DCMPLX( V1+V2,0.D0) 
PECO=RECO-v  CONST 
NBQ  =  0 

PRINT    4, <  DI( 1  ), 1=1, MAX) 
PRINT     19,NBQ,PEC0 
PRINT     14,NB2,VZ,DI ( MAX) 
DO    50  0    NP  =  NS,  NJ'4 
VQM    =    O.DO 
VI    =    C.CO 
V1J=O.DO 
V2    =    O.DO 
V2J=0.D0 
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DO    600     1=1, MAX 
J=MAX+1-I 

IF(DN(J,NP).LT.O.DO)     V1  =  VH-DN(J,NP) 

I F(DNJ( J, NP) .LT.O.DO)     V U = V 1J  +  DN J ( J , MP ) 

IF(DM(  J ,NP) .GE.O.DO)     V2=V2+DN(J,NP) 

IF( DM J(J,NP). GE.O.DO)     V2J=V2J+DNJ(J,NP) 

I F(DABS (DN(  J,NP) ) .  GT . DABS  I  V QM ) )     VQM=DN( J  ,  NP } 
600      CONTINUE 

P.EC(NP)  =  DCMPLX(  V1+V2,  V1J  +  V2J) 

REC(iVP)  =REC(NP)*CONST 

PRINT     19,N°,REC(  NP) 

PRINT  14,NP,VQM, DN(MAX,NP) 
50  0   CONTINUE 
725   CONTINUE 

VS=DCMPLX(O.DO,  l.DO) 

DO    800    NP=NS, NUM 
800       REC (MP ) =R EC( MP ) *VS**NP 

IF( IAVE.EQ.O )    GG    TO    960 

SUM=RECO 

IF(NS.ME.L)     SUM=O.ODO 

PRINT    19,NBQ,REC0 
C 

c 

C  CALCULATE    SUM    FOR    INTEGRATION    CHECK 


C 


DO    900    I=NS,NUM 

PQ=RFC( I) 

PRINT  19,I,REC(I) 
900   SUM=SUM+2.D0#PQ 

PRINT  18,SUM,SUMJ 

RETURN 
960   CONTINUE 

VS=DCMPLX(0.D0,-1.D0) 

PRINT    19,NPQ,REC0 

SUM=RECO 

IF(NS.ME.l)    SUM=O.DO 

SUMJ=O.DO 

DO    950    NP=MS,NUM 

VEC(NP)=REC(NP)*DCMPLX(DCOS(NP*-RANGL  )  ,DSIN(  NP*RANGL)  ) 

PRINT    19,NP,VEC(NP) 

PP=VEC(NP) 

SUM=SUM+2.D0*PP 

PP=VEC( NP)*VS 

SUMJ=SUMJ+PP 
950       CONTINUE 

PRINT    18, SUM, SUM  J 

RETURM 
285   KKK=1 

GO  TO  300 

1  F0RMAT110I5) 

2  FORMAT( 8E10.0) 

3  F0RMAT(2I10, 1P2E15.6) 

4  F0RMAT(1P8E15.5) 

5  FGRMAT(2I5) 

6  FORMAT ( 1P2E15. 6) 

12  FORMAT( 5X,G20.13) 

13  FCRMAT(5X,I4,6X,G20.13) 

14  FORMAT (5X,I4,6X,G20.13,6X,G20.13) 

15  FGRNAT( 2X , I  3 , 2X , 1 PE 15 . 6 ,2X , 1 PE1 5 .6 , 2X , I  3 , 2X , 1 PE15 .6  ,  2X 
1,  1PE15.6) 

16  F0RMAT(3(1P020.12 ) ) 

17  FORMAT(  5X,  ,V<AX    MOT     LARGE    ENOUGH     ZMAX  =  «  ,  Fl  5  .  8  ) 

18  FORMAT ( 2X, •SUM=« ,G20. 13, ' SUMJ=' ,G20. 13) 

19  F0RMAT(5X,I4,6X,G2  5.13, G20.13) 
END 
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FUNCTION  F  <R,TH,NN) 
C 
C 

C      FUNCTION    F  IS  AN  AUXILIARY  ROUTINE  cOR  SIGNAL  THAT 
C      COMPUTES  THAT  PART  OF  THE  INTEGRAND  NOT  RADIALLY 
C      SYMMETRIC 
C 


C 


IMPLICIT    REAL38     (A-H,0-Z) 

CCMPLFX»16    F 

D  I  MENS  ION    TH( 1) 

DIMENSION    A£E12^) , BBE(24) , A AO (24 ) , BBO( 24 ) 

COMMON     /SIPRM/    A,B,C,PX,QX 

DATA    NPR/24/ 

JJ=(~1)**NN 

AA=0.000 

BB=O.ODO 

IF(NN.GT.O)    GO    TO    5 

DO    1    J=  l^NPS 

X=R^DCOS(TH(J )) 

Y  =  R-VDSIN(  TH(  J)  ) 

DSX=DS I N(PX'X) 

DSY  =  DSI  N(  QXt  Y) 

DCX=DCOS(PX*X) 

DCY=DCOS(QX^Y) 

Pl  =  l.DO 

IF(X.NE.O.DO)     Pl=DSIN(A*X)/( A*X) 

IF(Y.NF.O.DO)     P1  =  P1-' DSIN(B*Y  )/  (  B#Y  ) 

AAE( J)=DCX*DCY*P1 

BBE(  J)  =  DSX---DSY-*-Pl 

AAO(  J)  =  -DCX    DSY'*P1 

BBO(  J)=-DSX**DCY*P1 
1  CONTINUE 

5  IF(JJ.LT.O)     GO    TO    10 

DO    20    J=1,NPR 

AA=AA+AAE(  J  ) '■  DCOS  (NN-*-TH(  J  )  ) 

B3=BB+BBE( J ) *DSI N ( NN?TH ( J ) ) 
20  CONTINUE 

F  =  DCMPLX( AA,  BB) 

IF(NN.cQ.O)     F=DCMPLX( AAtO.DO) 

RETURN 
10         DO    40    J=1,NPR 

AA=Afi+AAC( J) *DSIN(NN* TH  (  J  )  ) 

BB=BB  +  BBO( J)"  DCCS INN*TH(J ) ) 
40         CONTINUE 

F=DCMPLX(AA,BB) 

RETURN 

END 
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c 

c 
c 
c 
c 
r. 
c 
c 
c 
c 
c 
c 
c 


10 
99 


100 


SUBROUTINE    WAA1R1. R2, R,TH, W,N) 


SUBROUTINE    WAA    COMPUTES    THE     RADIAL 
ANGLES     AT    WHICH    THE     INTEGRAND    MUST 


NP    =     DEGREE    OF    GAUSSIAN    LA 
R    =    RADIAL    POINTS    TO    BE    USED 
TH    =    ANGLES    TO    BE     USED 
W    =    WEIGHTS    TO    BE    USED 
NFIRST     PREVFNTS     COMPUTING 


POINTS, WEIGHTS     AND 
BE    EVALUATED. 


G-L    POINTS    MORE    THAN    ONCE 


IMPLICIT    REALMS ( A-H.O-Z) 

DIMENSION    R(  1  ),*'(  1)  ,TH(1) 

REALMS    RC( 24), WC( 24) 

COMMON    /NFIRST/     NFIRST 

DATA    NP/24/ 

PI=DARCOS (-1.D0) 

PIN=PI /2.D0/NP 

IF(NFIRST,.NF.O)  GO     TO    99 

CALL    GL024(RC.^C0 .DO, 1 .DO) 

DO    10    I=1.NP 

TH( I )=( I-.5D0)*PIN 

NFIRST     =    1 

CONTINUE 

N    =    NP 

T    =    R2     -    Rl 

TW    =    T*-  (R2    +    Rl  ) 

F    =    TW^PI/NP 


DO 

RX  = 
RX  = 
RY  = 
RY  = 
R(  I  ) 
W(I  ) 


100 
=    Rl 


1=1 .NP 

'Rl 


RX- 

R2 

RY-.  R2 

=     DSORT(RX 

=     WC  {  I  M*  F 


+    RC(I)*(RY    -    RX)  ) 


CONTINUE 

RETURN 

END 
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SUBROUTINE    GL024 ( X , A, C , D) 
C 

c 

C  GAUSSIAN-LEGENQER     INTEGRATION    ROUTINE    (24    POINTS) 

C 

c 

DOUBLE     PRECISION     DMC, DPC 

DOUBLE    PRECISION    C,D,X(     24). A(     24),XX(     12),     AA(     12) 

DATA    XX(1),AM1).XX(2),AA(2).XX(3),AA(3),XX(4),AA(4),X 

1  XX(6)  ,AM  6)  ,XX(7) ,AA(7),XX(6),AA(6),XX(9),AA(9),XX(10 

2  XX(11 ) ,AA( 11) ,XX(  12) ,£A(  12)/ 

*.9951 67 2 1999  7 0  2 13  60 17  999  740 9 7 D-0, . 12  341 22  9 7 999  87 199546 
*. 97472  8 555  971 3 0949 819 83919  93 00-0,  .28  53138  8628933663181 
v.  938274 55 2002 73 275 8523 64900 17D-0, .442  7743  88174198  06163 
4  .8864 15 52 7 00440 103421 3 154 34 19 D-0,  .59298 5  849 154367  3  0 74o 
*.  8200  0 1  985  973  902  9  2  19  53  949  8  726D-0,.  73  3464  81 41  i  060  3  0  5  7  34 
*  .740124191578  5  5436424  38  28  10  310-0, .8619  01615319532  75917 
t. 648093 65 193 697 5 56 92 52495 7 8 69 D-0, .9  7613o5210411388  82  69 
*. 54 542 1471 38883 9 5 356 58375 61 72D-0,. 107444  270  1159  65  63473 
'.'.  43  3  79  3507626045  136  48  7  08  42  3190-0,.  115  50  5668C53  72  5601 35 
*. 3 15042 679696 163 37 43 66 7932 9 13 D-0, .121670472  9278  03  39120 
*. 19 1 1 18 86 747 36 16 309 15 8639 82 070- 0, .12  583  745  634632  8  2  9612 
4-.  64056  8  92  8626056  2  6  OS  5043  03  2  62  0-1,  .1279  3  8  19534675215697 

DMC    =     . 5D0^( D-Ci 

DPC    =     .5D0-MD+C) 

DO    2    1=1.12 

NI    =    25    -    I 

X(I  )    =    -DMC' XX ( I )     +    DPC 

X(NI)     =    DMC- XX ( I )     +    DPC 

A(I)     =     DMC»AA(I) 
2     MNI)     =     DMC-.  AA(  I  ) 

RETURN 

END 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 


c 
c 
c 
c 


c 
c 
c 
c 


SUBROUTINE    BESSL     ( AG* BF ,BZ , NUM) 


SURROUTINF    BFSSL    COMPUTES    BESSEL    FUNCTIONS    OF    THE 
RECEIVED    ARGUMENT 

AG    =    ARGUMENT    OF    BESSFL    FUNCTION 

BF    =    RETURN    VECTOR    OF    BESSEL    FUNCTIONS 

B/:    =    ZERO    ORDER     BESSEL    FUNCTION 

NUM    =     HIGHEST    ORDER    OF    BESSEL    FUNCTION     TO    BE    RETURNED 

NY. NX    =    ORDERS    CALCULATED    BY    ASYMPTOTIC     EXPANSION 


IMPLI 
DIMEN 
NZ  =  10 
KT  =  0 
IF(  AG 
I  F  (  AG 
X=AG 
NY=IF 
I  F  ( NY 
I  F  (  NY 
IF  (NY 
I  F  ( NY 
IF(NY 
IF(NY 
I  F  (  NY 
IF  (NY 
I  F  ( NY 
IF(NY 
IF(NY 
I  F(  IF 
IF(NY 
NY=NF 
NX  =  NY 


CIT  REAL^8  (A-H.O-Z) 
SIGN  BF3(  1000)  ,BF(1) 
00 

.LT.O)     KT=1 

.LT.O.DO)     AG=DABS(AG) 


IX 
.L 
.E 
.F 
.G 
.G 

•  G 
.G 
.G 
.G 

•  G 
.G 
IX 
.L 

-1 


(SNGL( 
E.l)  N 
0.2)  N 
0  .  3  )  N  F 
E.4.AN 
F.6.AN 
E.10.A 
E.15.A 
E.29.A 
E.80.A 
E.150. 
E.330) 
(FLOAT 
T.l)     N 


AG)  ) 
F=16 
F  =  20 
=  2<+ 

D.NY. 
D.NY. 
ND.NY 
ND.NY 
NO.  NY 
ND.NY 
AND.N 
MF=I 
(NF)/ 
F  =  8 


LT.6)    NF  =  6=?NY 
LT.IO)        NF=5*NY 
.LT.15)     NF=4^NY 
.LT.29)     NF=3*NY 
.LT.80)     NF=2*NY 
.LT.150)     NF=IFIX(1. 
Y.LT.330)     NF=IFIX(1 
FIX(  1.25vFL0AT(NY) ) 
2.) .LT.IFIX(FLGAT(NF+ 


5i FLOAT(NY) ) 
.4'-v  FLOAT  (NY)  ) 


l)/2.))     \F=NF+1 


SET  ORDERS  >  NY  TO  ZERO 

DO    1    I=NY.NZ 
BF3(I  )  =  0.D0 

CONTINUE 


COMPUTE    ASYMPT3TIC    EXPANSION 


DO  2 

z=x/ 

AMU  = 
A  = 
Al=l 
DEL 
XI  = 
CALL 
Dl  = 
D2  = 
BO  = 
IDO-rA 
FI  = 
BF3( 
CONT 
NY2  = 


MU=NX 

DFLOAT 

DFLOAT 

DSORT( 

.5D0^( 

=    DEXP 

DEXP( 

TABLE 

DEXP( 

DEXP( 

-5. DO 

<y3)-l 

DEXP( 
MU) =  F-I 
INUE 
NY-2 


.NY 
(MU) 
(MU) 
l.DO 
DLGG 
(2.D 
2. DO 
(XI 
l.DO 
5. DO 
/(43 
.ODO 
0.25 
M  AI 


(  (  1. 

0/3. 

/3.D 
.AIV 
/3.D 
/3.D 

.DO?- 

/( a. 

DO  :<D 
\//Dl 


2) 

DO 
DO 
0* 
♦  A 
0* 
0* 
DE 
OD 
LO 
+  A 


+A)/Z)-A) 
tDLOG(Al)) 

DLOG(AMU)+DLGG(DEL) ) 

IPV) 

DLOG(AMU) ) 

DLOG(AMU) ) 

L»-*2)  +  (  1.D0/DSQRT(DEL  )  )-H5 

0*  A )  ) 

G(  4.  DO*- DEL  )-0.5D0>DL0G(A)  ) 

IPV-B0/D2) 


D0/(24. 


USE    RECURSIVE    EQUATION    TO    GENERATE    ORDERS    <    NX 

DO    4    I-1.NY2 

N  =  i\Y2-l+l 

BF3(N)     =     (2.D0-DFL0AT(N  +  1 )  /  X  )  .-3  F3  I  N  +  l )-BF3(N+2) 
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c 
c 

c 
c 


4     CONTINUE 

SUM  =  2.D0-rBF3(  1  )/X-BF3(2) 

DG  76  I =2.NY2,2 

SUM  =  SUM*2.D0.*3F3(I) 

76    CONTINUE 


C 

c 
c 
c 


COMPUTE  NORMALIZING  CONSTANT 

SNORM  =  1. DO/SUM 

DO  77  1=1, NUM 

BFl I) =SNORM^BF3 ( I) 


ODD    ORDER    BESSEL    FUNCTIONS     ARE    ODD 

IF(KT  .EO.l)     BF{  I )  =  BF( I )* ( - l.DO) $+ I 
7  7  CONTINUE 

BZ    =    ( 2.DO*SF3< 1)/X-BF3(2) J* SNORM 

RETURN 

END 
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c 
c 
c 
c 
c 
c 


2 

4 

6 
8 


12 

13 

19 


14 
30 


SUBROUTINE     TABLE ( X , AI V. AI PV) 


SUBROUTINE    TABLE     IS    AN    AUXILIARY    ROUTINE    FOR    BESSL 
THAT    COMPUTES    AIRY    FUNCTIONS 


IMPLIC 

DIMENS 

DATA  A 
1-55872 
2.55391 
3. 54958 

DATA  A 
1. 57192 
2. 5788  7 
3.58  52  3 

DATA  A 
1. 54084 
2.53060 

DATA  A 
1.59837 
2.61427 

L  =  0 

Y  =  DE 
W  =  DE 
L=l 
IF(Y. 
I  F(  Y. 
L  =  2 
Wl  = 
I W1=I F 
W2=W1+ 
IW2=IF 
I  Z  =  I  W 1 
IF( IW1 
IY=IZ- 
IF( IY. 
GO    TO 

Y  =  10 
I Y=IFI 
IF(  IY- 
I Z=IY+ 
GO  T08 
AIV=AI 
AIPV=A 
GO  TO 
IZ=IY 

I Y=IYt 
DP=Y-F 
IF(L.E 
I  F  (  L  .  E 
AIV  = 
AIPV  = 
GO  TO 
AIV  = 
AIPV  ^ 
GO  TO 
DP  =  Y 
AI  V  = 
AIPV  = 
AIV  = 
AIPV  = 
AIPV  = 
RETURN 
PRINT 
FORMAT 
RFTURN 
END 


IT  REAL 
ION    AH 

1  /4*0.0 
4D0,  4!   0 

2  DO. 41 0 
4D0.4r 0 
I  P/4-.  0. 
7  DO .4*0 
800.  4*0 
5  CO  ,  4~  0 
B 1/4-0. 
400. .53 
1  DO.  .52 
ttIP/4*0 
200, .60 
5  DO  .  .  6  1 

XPCDLOG 
XP(DL0G 

T .1 .500 
T. 0.499 

0  .DO    •*■ 
I X( SNGL 
1. 

IX (SNGL 
*10 

.LT-IW2 
5 

EO.O)  G 
8 

.DO  •*  Y 
X( SNGLt 
Y )  2.4, 
1 

(  IY  ) 

1  PC  IY) 
19 


*8 
50) 

DO. 
.OD 
.OD 
.OD 
OCO 
.OD 
.OD 
.00 
ODO 
861 
878 
.00 
178 
715 

(  1 
(2 

)  G 
999 

Y 
(Wl 

(W2 

)  I 
0  T 


Y)  ) 
6 


( A-H.O- 
, AIPC50 
.562280 


5570 
.5524 
.  5482 
56687 
.5743 
.  5810 
.5872 
,  .54823 
SDO, . 53 
3  DO. .52 
0..5872 
2D0.  .  60 
6D0, .61 


0, 
0. 
0. 

«    . 

0. 
0. 
0, 


Z) 

),  ABK  1 
DO. 4*0. 
58D0.4v 
21D0.4* 
3000/ 
300.  4-5-0 
20  00,4* 
5600,4* 


DO/ 


GDO.  .54 
6489D0, 
702700/ 
45D0  ,  .5 
5068D0. 
9954D0/ 


5)  , ABIPC  15) 

ODO.  .560  46  2DO,4'-.-0.ODO, 
0. ODO.. 5  5  545600,4^0. ODO, 
0.0  DO.  .55  0  98000.  4*'  0.  ODO, 

.ODO.. 56944 800,4*0.  ODO, 
0.0  CO,  . 576635 DO, 4 *0. ODO, 
0. ODO, .583 174 DO, 4* 0.0 DO, 

5636D0, .543180D0, 

. 534448D0, . 532488D0, 

91120D0, .594823  DO, 
. 608239D0, .611305D0, 


5D0)-1.5D0*DL3G( X) ) 
D0/3.D0J+1 .5D0*DL0G(X) i 

0  TO  14 

DO)  GO  TO  1 


)  ) 
)  ) 

Z=IZ+5 

0  13 


1 

LOA 
0.2 
0.1 
(  1. 

(1 
19 
(  I. 

(1 
19 
/5. 
(  1. 

(  1 

oex 

Or 
A I 

30 
(  5X 


T(  IY 
)  DP 
)   G 

DO-D 
.00- 

DO-D 
,.  DO- 
DO 

DO-D 
.00- 
P(-0 
XP(0 

pv  - 


) 

=  DP 
0  T 
Pit 
DP) 

P)«- 

DP) 


P)  * 
DP) 
LOG 
.25 
(- 


/5.D0 

0  12 

AI ( IY)+DP>AI(IZ) 

*AIP( IY)+DP~AIP( IZ) 

ABK  IY)+DP*ABI<  IZ) 
*ABIP(  IYi+DP«ABIP(  IZ) 


0.564190D0  +  DP  -Alt  IZ) 
*0. 564190D0+DPvA IP( IZ) 
(2.D0)-0.25D0-DLCG(X)-W  +  DL0G( A I  V)  ) 
DO*DLGG(X)-W+DLOG(AIPV) ) 
.5D0) 


•HAVE  EXCEEDED  TABLE  VALUES') 
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SUBROUT INF    SUMCK 
C 

c 

C  SUBROUTINE  SUMC<  COMPUTES  THE  ANALYTICAL  EXPRESSION 

C  FOR  THE  INFINITE  SUM  OF  COEFFICIENTS 

C 

C  SISUM  =  INFINITE  SUM  OF  POINT  TARGET  (IAVE=1)  OR 

C  AVERAGED  TARGET  (IAVE=0)  COEFFICIENTS 

C  BKSU^I  =  INFINITE  SUM  OF  BACKGR3UND  COEFFICIENTS 

C 

c 
c 

IMPLICIT    REAL*3     (A-H.C-Z) 

COMMON     /DTPRM/    RHO , SI G , WI » H , AL , BE . RX t RY , XO , YO 
COMMOM     /AVPRM/     .MUM,  I  AVE,  PO 
COMMON    /ANGL/     RANGL 
W  =  WI 

CONST=0.2500D0 
GAMA=DSORT(SIG^  *2+PO**2) 
IF(  IAVE. EO.O)     GAMA  =  SI G 
DIV=2.D0*DSQRT(2.D0)*GAMA 
IF{ IAVE. EO.O)     GO    TO    10 
C 

c 

C      COMPUTE  SUM  OF  AVERAGED  COEFFICIENTS 
C 

A=DERF(  (2..D0MX0  +  RHO)  +  W)/DI  V) 

B=DERF( (2^DO*(X0+RHO)-WJ/DIVJ 

C=DERF<  <2.0D0<M  YO)+H)/DIV) 

D=DERF( (2.0D0*(Y0)-H) /DIV) 

4  SI  SUM  =  CONST'f  i  A-  B  )  >  (C-D) 
R1=R(  WtAL,SIGJ 
R2=R(H, BE.SIG) 

C 

C 

C      COMPUTE  SUM  OF  3ACKGROUN0  COEFFICIENTS 

C 

BKSUM=4.D0*AL*BE*R1*R2/(W*H)**2 
PRINT  5.BKSUM.SISUM 
RETURN 
10    CONTINUE 

C 

C 

C      COMPUTE  SUM  OF  POINT  TARGET  COEFFICIENTS 

C 

A  =  DERF(  (2-0 DO* ( XO+RHO-DCOS ( RANGL  J-RXJ+W  J /DIV) 
B  =  DERF(  (2.0D0*(  X0+ RHO*  DCGS  (  RANGL  ) -RX  ) -W  )/  CIV) 
C=  D ER  F  (  (  2  .  DOC  (  Y  0  +  RHO"'  D  S  I  N  (  R  A\  3  L  )  -R  Y )  +  H )  /  D  I  V  ) 
D =DE  RF  (  ( 2  .  DO*-  (  Y 0  +  RHO*  DS  I N  {  R  AN GL  )  -R  Y  )  - H  )  /  D  I  V  ) 
GO    TO    4 

5  FORMAT ( 2X.«BKSJM=',D15.8, 10X,  • SISUM=« .D15.6) 
FND 
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c 
c 
c 
c 
c 


FUNCTION    R(U.Z.SIG) 


FUNCTION    R    IS    AN    AUXILIARY     ROUTINE    FOR    BKSUM 

IMPLICIT    REAL*8     (A-H.O-Z) 

PI=DARCOS(-1.DO) 

R=O.DO 

Fl=-< (U/(2.D0*SIG)  )**2  ) 

E1=DEXP(F 1) 

E2«(Z*S  IG)-:-'>2 

E2=DEXP(E2) 

E31=Z*U 

E3=DEXP (E31 ) 

E4=DEXP(-E31 ) 
X1  =  U/ ( 2  .DOvSIGJ 
Yl=l.DO-DERF(Xl ) 
X2=Z£SIG+U/<2.D0SSIG) 
Y2=I .D0-DERF(X2) 
X3^Z*-SIG-U/(2.D0->S  IG) 
Y3=1.D0-DERF( X3) 
X4=Z*SIG 
Y4=l.DO- DERF(X4) 

R=(U/(Z**2)  )<■  (  1.00-Yl-(  SI  3s:  2.DO/(DSORT(  PI  )*U)  )*(1.D0-E 
R=R+(E2/(2  .D0*Z**3J  I*  (  E3*Y  2  +  E4*-  Y3-2.D0*  Y4) 
RETURN 
END 
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